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A THEORETICAL STUDY OF MIXING DOWNSTREAM OF TRANSVERSE 


INJECTION INTO A SUPERSONIC BOUNDARY LAYER 
A. J. Baker and S. W. Zelazny 
Bell Aerospace Company 


SUMMARY 


This report presents the results of a theoretical and analytical study of mixing downstream of 
transverse hydrogen injection, from single and multiple orifices, into a Mach 4 air boundary layer flow 
over a flat plate. Numerical solutions to the governing three-dimensional, elliptic boundary layer equa- 
tions were obtained using a general purpose computer program (COMOC) founded upon a finite ele- 
ment solution algorithm. The concepts of two-dimensional turbulent mixing length and mass defect 
theories were extended to establish a prototype three-dimensional turbulent transport model. 

Excellent agreement between the COMOC-computed flow field and experimental data for a 
jet/freestream dynamic pressure ratio of unity was obtained in the centerplane region of the single jet 
configuration. Poorer agreement off centerplane may be indicative of the inadequacy of the extrapo- 
lated two-dimensional turbulence model that neglects flow field three-dimensionality. A considerable 
improvement in off-centerplane computational agreement occurred for a multi-jet configuration, using 
the same turbulent transport model. 

Some of the computational disagreement with measurements may well reflect initial condition 
variability. One alternative to requiring detailed initial data is to employ the “virtual source” con- 
cept, wherein the complex transverse injection phenomena is computationally replaced by a hydrogen 
jet imbedded within the air boundary layer, the distinguising features of which are solely a function of 
freestream and injectant parameters. A theoretical model for establishing initial conditions for a virtual 
source was derived. Very good computational agreement with experimental data was obtained for the 
multiple injector geometry for the three dynamic pressure ratios investigated (q r = 0.5, 1 .0, 1 .5), at 
stations up to 60 (injector orifice) diameters downstream. Single jet comparisons could be obtained 
using the virtual source concept, if desired by a simple change in lateral boundary condition specifi- 
cations. 


The COMOC-computed solutions for the three-dimensional flow field provide considerably de- 
tailed information, much more than is typically available by experimental means. Important design 
guidance for the engineer may be obtained by integrating these data over the problem solution domain 
to generate scalar correlating parameters. Since the prediction of reacting flows is an ultimate goal, a 
simple integral parameter, defined as the percentage of hydrogen that could stoichiometrically react in 
a given concentration and velocity distribution within the air boundary layer was computed. Shown in 
figure 1 are the COMOC-computed mixing efficiencies, 17 , for the “virtual source” modeling of multi-jet 
injection, as a function of distance downstream from the injection orifice. Mixing efficiency is computed 
to be noticeably poorer for q r = 1.5 than the other two dynamic pressure ratio conditions. 

COMOC can be readily expanded to predict either equilibrium or finite rate chemistry for the 
three-dimensional mixing problem. On a relative basis, however, if reaction is fast, hence diffusion con- 




Figure 1. COMOC Computed Mixing Efficiency for Virtual Source, Multi-Jet 
Configuration for Three Dynamic Pressure Ratios 

trolled, the computation of mixing efficiency for a cold flow configuration may provide important de- 
sign guidance, by capturing the detailed three-dimensional computed flow field description within a 
single scalar parameter. 


INTRODUCTION 


The hydrogen fueled scramjet engine is a prominent candidate for propulsion of the next gen- 
eration of hypersonic cruise vehicles, see for example Becker and Kirkham (ref. 1) and Bushnell (ref. 2). 

An airframe-integrated underbody engine configuration has been suggested (ref. 3), and design con- . 
siderations are discussed by Henry and Anderson (ref. 4). Over the years many alternative designs have 
been proposed by the Air Force, Navy, and NASA as well. They all enjoy a certain commonality in 
that fuel introduction arrangements typically consist of rows of circular, choked flow fuel injector 
orifices mounted flush or normal to the combustor wall or in fins spanning the combustor inlet. The 
various proposed component designs have largely emerged from experimentation wherein empirical 
relations have established a preliminary configuration for a starting point. Early attempts to estab- 
lish a theoretical basis for this developing technology have been confined by lack of detailed data 
over the complete three-dimensional flow field. Since the pattern of fuel injection, hence three-dimen- 
sional mixing, will exert significant influence on combustor performance, an in-depth analytical analysis of 
the complex three-dimensional flow fields is required to lend critical support to design technology. 

The results of comprehensive experimentation on prototype injector subsystems are now 
becoming available. For example, Rogers has performed detailed measurements over the entire 
three-dimensional mixing region downstream of hydrogen injection from a single discrete orifice 
(ref. 5), and multiple, laterally-disposed orifices (ref. 6), into a Mach 4 airstream over a flat plate. 

Similar measurements for a Mach 2.7 airstream, including Schlieren photographs of the near injec- 
tion region, are reported by Wagner, Cameron and Billig (ref. 7). These data are not only useful for 
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over a flat plate. Numerical solutions to the governing three-dimensional, elliptic boundary layer equa- 
tions were obtained using a general-purpose computer program (COMOC) founded upon a finite ele- 
ment solution algorithm. The concepts of two-dimensional turbulent mixing length and mass defect 
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Excellent agreement between thmCOMOC-computed flow field and experimental data for a 
jet/freestream dynamic pressure ratio of unity was obtained in the centerplane region of the single jet 
configuration. Poorer agreement off centerplane may be indicative of the inadequacy of the extrapo- 
lated two-dimensional turbulence model that qeglects flow field Three-dimensionality. A considerable 
improvement in off-centerplane computational\greement occurred for a multi-jet configuration, using 
the same turbulent transport model. ' 


Some of the computational disagreement with measurements may well reflect initial condition 
variability. One alternative to requiring detailed initial data is to employ the “virtual source” con- 
cept, wherein the complex transverse injection phenomena is computationally replaced by a hydrogen 
jet imbedded within the air boundary layer, the'tiistinguiting features of which are solely a function of 
freestream and injectant parameters. A theoretical model for establishing initial conditions for a virtual 
source was derived. Very good computational agreement with experimental data was obtained for the 
multiple injector geometry for the three d/namic pressure radios investigated (q r = 0.5, 1 .0, 1 .5), at 
stations up to 60 (injector orifice) diameters downstream. Single jet comparisons could be obtained 
using the virtual source concept, if desired by a simple change in\ateral boundary condition specifi- 
cations. 

The COMOC-computed solutions for the three-dimensional flW field provide considerably de- 
tailed information, much more than is typically available by experimental means. Important design 
guidance for the engineer may beobtained by integrating these data over^the problem solution domain 
to generate scalar correlating parameters. Since the prediction of reactingvflows is an ultimate goal, a 
simple integral parameter, defined as the percentage of hydrogen that could\stoichiometrically react in 
a given concentration and velocity distribution within the air boundary layenwas computed. Shown in 
figure 1 are the COMOC-computed mixing efficiencies, 17, for the “virtual source” modeling of multi-jet 
injection, as a function ofi^listance downstream from the injection orifice. Mixing efficiency is computed 
to be noticeably poorenTor q = 1.5 than the other two dynamic pressure ratio conditions. 

/ v 

COMOC can^be readily expanded to predict either equilibrium or finite rate chemistry for the 

three-dimensional mixing problem. On a relative basis, however, if reaction is fast, hence diffusion con- 

/ \ 


determining the mixing resulting from various given injector configurations but, of equal importance, 
they provide the required accuracy comparison for theoretical analysis of the resultant turbulent flow 
fields. Concomitant with the establishment of appropriate models for turbulent mass and momen- 
tum transport coefficient distributions in three-dimensional flow fields will be the emergence of the 
capability to perform analytical studies for both detailed flow description and design optimization. 
The success of this approach depends upon establishment of a proven capability to model the three- 
dimensional viscous flow field equations numerically, and to generate computational solutions 
for practical geometrical configurations. 

This numerical analysis capability is now operational within the COMOC general purpose 
computer program system. The theoretical foundation for COMOC is a finite element solution al- 
gorithm for a generalized initial-valued, quasi-linear elliptic boundary value partial differential equa- 
tion, systems of which typically describe the mechanics and thermodynamics of continuum mech- 
anical phenomena. Employing an extrapolation of two-dimensional mixing length/mass defect tur- 
bulent transport concepts, as a prototype three-dimensional turbulent transport model, an analy- 
tical study of mixing downstream of transverse hydrogen injection into a Mach 4, three-dimen- 
sional air boundary layer flow over a flat plate has been performed using COMOC. Detailed num- 
erical prediction of the flow field is presented, in a solution domain spanning 30 to 1 20 injector 
diameters downstream, for the single jet geometry of Rogers (ref. 5), for a jet to freestream dynamic 
pressure ratio (q r ) of unity, an injection orifice diameter (D) equal to 0.102 cm (0.040 in.), stagnation 
freestream pressure of 13.6 atm (1 atm = 101.325 kN/m 2 ) and a stagnation temperature of 300°K. 
These conditions correspond to a Reynolds number per meter of 6. 19 x 10 7 ; the boundary layer 
thickness at the injector longitudinal station was equal to 2.7 injector diameters. The results of com- 
putations for a multi-jet configuration, with injector separation equal to 12.5 injector diameters, are 
also presented, for the freestream conditions similar to the single-jet case. 

Considerable experimental data management is required to establish initial conditions for 
numerical three-dimensional solutions. A candidate alternative method is computational replace- 
ment of the complex injection and turning phenomena with a “virtual source” of injectant imbedded 
within the air boundary layer flow. A theory is presented to establish the initial condition specifi- 
cation of a virtual source configuration based solely upon freestream and injectant parameters. Em- 
ploying the same extrapolated turbulence model as before, analytical results are presented for pro- 
pagation of the virtual source of hydrogen into the downstream mixing region, for a solution do- 
main spanning 0 to 60 injector diameters downstream. The computational configuration is identical 
to the multi-jet geometry of Rogers, and results are established for the three dynamic pressure ratios 
of q r = 0.5, 1.0, and 1.5. The agreement between computed and experimentally measured data is 
described in terms of decay of the maximum hydrogen concentration, trajectory of the maximum 
concentration, and lateral spreading of the diffusion pattern. 

For all computational solutions, COMOC was adapted to compute the integral mixing effi- 
ciency parameter, 17, of ref. 6. The potential usefulness of this scalar correlation parameter as a 
design optimization tool on a relative basis is discussed. Combu&on data are required to establish 
an absolute comparison basis. Sample results are discussed illustrating how the COMOC predictive 
analysis capability can be expanded to compute the actual reacting flows of practical interest in 
combustor design and optimization. 
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I 


THE DIFFERENTIAL EQUATION SYSTEM GOVERNING THREE- 
DIMENSIONAL MIXING FLOW FIELDS 

The description of a state point in viscous fluid mechanics is contained within the solution 
of the system of coupled, nonlinear, second-order partial differential equations enforcing local con- 
servation of species mass, total mass, linear momentum, and energy. Closure of this equation system 
requires identification of constitutive laws. For laminar flows, transport properties such as viscosity 
and thermal conductivity are describable in terms of molecular behavior and therefore dependent 
solely upon the material present. In turbulent flows, the time-averaged Navier-Stokes equations 
appear identical to the laminar flow equations, after identification of turbulent or “eddy” transport 
coefficients which are dependent upon the kinematics of the flow field and not the molecular prop- 
erties of the fluid. In this development, a generalized transport property description is assumed, 
which may be selectively laminar or turbulent as required. 

An illustration of the subject three-dimensional flow field is shown in figure 2. The pre- 
dominant direction of flow is aligned with the x axis. Although the essential character of the flow 
is boundary layer, the common three-dimensional boundary layer equations, (ref. 8) are an insuffi- 
cient description, since some dependent variable gradients in the lateral direction (z axis) may not be 
uniformly negligible (in comparison to those parallel to the y axis). The governing equation system 
is thus the elliptic boundary layer form of the parabolic Navier-Strokes equations, which retain all 
derivatives in the lateral flow direction. The three-dimensional velocity vector U has scalar com- 
ponents in the rectangular Cartesian basis, figure 2 as 

— » A A A 

U = ui+vj+wk (1) 

It is useful to identify a two-dimensional vector differential operator V 2 , 

V 2 Sj( ),y + k( ), z (2) 

and the comma notation signifies partial differentiation. The three-dimensional, elliptic boundary 
layer equations for a multispecies and compressible fluid can then be written as, 

Global Continuity: 

(3) 


(4) 


puu, x =V 2 '(mV 2 u) -pU«V 2 u -p, x (5) 

Lateral Momentum: 

puw, x = V 2 .(pVjw) - pTT* V 2 w -p, z (6) 


0=V 2 .(pU) +(pu), x 

i^ Species Continuity: 

puY| x =V 2 .(^ e v 2 Y i )-pU-V 2 Y'+S' 

A p r 

Longitudinal Momentum: 
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determining\he mixing resulting from various given injector configurations but, of equal importance, 
they provide \he required accuracy comparison for theoretical analysis of the resultant" turbulent flow 
fields. Concomitant with the establishment of appropriate models for turbulent mass and momen- 
tum transport coefficient distributions in three-dimensional flow fields will be thenimergence of the 
capability to perform analytical studies for both detailed flow description and design optimization. 
The success of this approach depends upon establishment of a proven capability to model the three- 
dimensional viscous flow field equations numerically, and to generate computational solutions 
for practical geometrical configurations. 

\ / 

This numerical analysis capability is now operational within the COMOC general purpose 

computer program system. Vhe theoretical foundation for COMOC is a finite element solution al- 


gorithm for a generalized initial-valued, quasi-linear elliptic boundary value partial differential equa- 
tion, systems of which typically describe the mechanics and thermodynamics of continuum mech- 
anical phenomena. Employing an extrapolation of two-dimensional mixing length/mass defect tur- 
bulent transport concepts, as a prototype three-dimensional turbulent transport model, an analy- 
tical study of mixing downstream of transverse hydrogen injection into a Mach 4, three-dimen- 
sional air boundary layer flow over a\flat plate has been performed using COMOC. Detailed num- 
erical prediction of the flow field is presented, in a solution domain spanning 30 to 120 injector 
diameters downstream, for the single jet geometry of Rogers (ref. 5), for a jet to freestream dynamic 
pressure ratio (q r ) of unity, an injection orifice diameter (D) equal to 0.102 cm (0.040 in.), stagnation 
freestream pressure of 13.6 atm (1 atm = Ml. 325 kN/m 2 ) and a stagnation temperature of 300°K. 
These conditions correspond to a ReynoldsVumb/r per meter of 6. 19 x 10 7 ; the boundary layer 
thickness at the injector longitudinal stationVa^equal to 2.7 injector diameters. The results of com- 
putations for a multi-jet configuration, with injector separation equal to 12.5 injector diameters, are 
also presented, for the freestream conditions,similar to the single-jet case. 

/ \ 

Considerable experimental data management is required to establish initial conditions for 
numerical three-dimensional solutions. candidate, alternative method is computational replace- 
ment of the complex injection and turning phenomena with a “virtual source” of injectant imbedded 
within the air boundary layer flow. A^i theory is presented to establish the initial condition specifi- 
cation of a virtual source configuration based solely upqn freestream and injectant parameters. Em- 
ploying the same extrapolated turbtilence model as before, analytical results are presented for pro- 
pagation of the virtual source of Jiydrogen into the downstream mixing region, for a solution do- 
main spanning 0 to 60 injector diameters downstream. TlAcomputational configuration is identical 
to the multi-jet geometry of Rogers, and results are established for the three dynamic pressure ratios 
of q r = 0.5, 1.0, and 1.5. The agreement between computed and experimentally measured data is 
described in terms of decaff the maximum hydrogen concentration, trajectory of the maximum 
concentration, and lateral/spreading of the diffusion pattern. \ 
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For all computational solutions, COMOC was adapted to compute the integral mixing effi- 
ciency parameter, 77, of ref. 6. The potential usefulness of this scalar^orrelation parameter as a 
design optimization ^t'ool on a relative basis is discussed. Combu3ion d*ata are required to establish 
an absolute comparison basis. Sample results are discussed illustrating how the COMOC predictive 
analysis capability' can be expanded to compute the actual reacting flows^f practical interest in 
combustor design and optimization. \ 
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Figure 2. Three-Dimensional Flow Field for Mixing Downstream of 
Transverse Injection from Discrete Orifices 

Energy: 

+ w 2 )] 

Y*] (7) 

The variables appearing in equations (3) through (7) have their usual fluid mechanic interpretation: 
Y 1 is the mass fraction of the ith species, Pr is the Prandtl Number, and Le is the Lewis Number 
specified for essentially binary diffusion. The stagnation enthalpy H is defined as, 

H = Z h' Y* + 1/2 (u 2 +w 2 ) 
i 

and the static enthalpy of the ith species is expressed in terms of temperature as, 

T i 

h‘=f CpdT 

Each species is assumed to behave as a perfect gas. From Dalton’s law, the equation of state takes 
the form 

Y 1 

p = pRT 2 — (10) 

i W 1 

where R is the universal gas constant, and W 1 is the molecular weight of the ith species. 

Solution is required for the equation system (3) through (10), with appropriate specification 
of both boundary and initial conditions, as well as the pressure distribution, p (x,z), in the external 
inviscid flow field. The form selected for writing these equations purposefully illustrates the 
mathematical uniformity that pervades the system which is of considerable significance in establish- 
ing the (finite element) numerical solution algorithm. 




( 8 ) 

(9) 


puH=V 2 -(~V 2 H) -pU-V 2 H 
Pr 

-V 2 .f(l - Pr) £ V 2 (u 2 
2Pr 

-V 2 .[(l - Le) Eh i V 2 



5 




FINITE ELEMENT SOLUTION ALGORITHM 


Each of the partial differential equations (4) through (7), requiring solution for the three- 
dimensional mixing problem, is a special case of the general initial-valued, quasi-linear elliptic 
boundary value problem of mathematical physics. The global continuity equation (3) is purely 
initial value, and requires special consideration as detailed after the general development. For the 
numerical algorithm development, full advantage is taken of the mathematical uniformity pervading 
equations (4) through (7) by noting that each is a special case of the general equation, L (q) = 0, or 
specifically, 

puq, x = V 2 -(kV 2 q)-pU- V 2 q + f (11) 

In equation (1 1 ), q is a generalized dependent variable, k is the generalized diffusional transport co- 
efficient, and f is any forcing function that may or may not be explicitly independent of q. Table 1 
lists the functions k and f appropriate to q identified with each dependent variable. 


TABLE 1 

COEFFICIENTS IN GENERAL EQUATION 


Equation 

Number 

q 

k 

f 

(4) 

y‘ 

pLe 

Pr 

s‘ 

(5) 

U 

P 

-P’x 

(6) 

W 


'P’Z 

(7) 

H 

P 

Pr 

-V 2 .[(l - Pr) iL V 2 (u 2 + w 2 )] 
2Pr 




-V 2 .[(l - Le)JLS h' V 2 Y'] 


The initial conditions for solution of equations (4) through (7) comprise a specified distri- 
bution for each q at the initial longitudinal station, x=x Q , of the form 

q(x 0 ,y,z) = Q 0 (y,z) (12) 

All boundary conditions of practical importance, for the elliptic solution domain, can be written 
in the form 

a, (x)q (y,z) + a 2 (x)v 2 Q G',z)-n = a 3 (x) (13) 

The superscript bar constrains the independent variable to lie on the closure, 3R, of the elliptic 
domain R, and n is the outward pointing unit vector everywhere normal to the domain closure. 

Table 2 lists representative values of the aj, equation (13), required to enforce sample boundary con- 
dition specifications for select dependent variables. 


6 





Flat Plate 


0 Injection Orifice 



Figure 2. 'Three-Dimensional Flow Field for Mixing Downstream of 
Transverse Injection from Discrete Orifices 


Energy: 


puH, x = V 2 • (f V 2 H)\ -pU-V 2 H 
Pr \ 

-V 2 .[(l -Pr\^V 2 (u 2 + w 2 )] 
\ r / 
-V a .[(l-U)|SWv 2 Y*/ 



The variables appearing in equations (3) through (7J) have their usual fluid mechanic interpretation: 
Y 1 is the mass fraction of the ith species, Pr is tlie^Prandtl Number, and Le is the Lewis Number 
specified for essentially binary diffusion. The stagnation enthalpy H is defined as, 


H = 2 h' Y 1 + 1/2 (u 2 +w 2 ) 
i 


and the static enthalpy of the ith species is^expressed ii\terms of temperature as, 


]i s i v 



Each species is assumed to behave as a perfect gas. From DaltoVs law, the equation of state takes 
the form / \ 

p = pRT 2 — / \ 

i W 1 / \ 

where R is the universal gas constant, and W 1 is the molecular weight of\he ith species. 

7 \ 

Solution is required for the equation system (3) through (10), with\appropriate specification 
of both boundary and^nitial conditions, as well as the pressure distribution ,\ (x,z), in the external 
inviscid flow field. The form selected for writing these equations purposefullyullustrates the 
mathematical uniformity that pervades the system which is of considerable significance in establish- 
ing the (finite element) numerical solution algorithm. \ 
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TABLE 2 

COEFFICIENTS IN GENERAL BOUNDARY 
CONDITION STATEMENT 


Boundary Condition 

Coefficients 

ai 

a 2 

a 3 

No-Slip at Wall 

1 

0 

0 

Slip at Wall 

t 

-1 

0 

Mass Injection 

0 

1 

t 

Adiabatic Wall 

0 

1 

0 

Specified Heat Flux 

0 

1 

t 

Temperature Dependent Flux 

t . 

t 


Symmetry Condition 

0 

1 

H 


t Non-zero values, user specified to enforce desired 
conditions. 


Solution of the differential equation system, equations (4) through (7) is thus reduced to 
establishment of a numerical solution algorithm for equations (11) through (13). The finite element 
method is suggested, as has been presented by Baker (ref. 9, 10). Briefly, solution is sought to 
equation (11) using Galerkin criteria within the Method of Weighted Residuals (MWR), see Finlay- 
son and Scriven (ref. 1 1 ) by approximation of a (each) dependent variable within the solution 
domain, R, by a series expansion of the form. 

q* (x) = 2 C k (x)0 k (y,z) (14) 

k = 1 

In equation (14), the C k (x) are unknown expansion coefficients modifying members, 0 k (y,z), of 
a set of functions that are complete in the space of the elliptic domain. Certain of the expansion 
coefficients are evaluated by requiring that the approximate solution, equation (14), satisfies ap- 
propriate boundary conditions, equation (13). The remaining unknown coefficients are then de- 
termined by setting to zero the differential equation residuals, formed by substitution of q* into 
equation (1 1), multiplication by a set of weighting functions, W k , typically identical to the spatial 
approximation functions, d> k , and integrating over the elliptic domain R. This produces the N 
ordinary differential equations. 

£w k (x) L (q*) dr = 0 k=l,2,...,N (15) 

where N corresponds to the total number of remaining unknown coefficients in the definition of 
the approximation function, equation (14), and L(q*) is equation (11) written on the approximate 
solution. 

When non-vanishing gradient boundary conditions exist, a generalization to MWR relaxes 
the constraints formed by the boundary condition statement. Form the N boundary residuals, de- 
fined from equation (13) as 










Looking to the variational calculus for guidance, multiply equation (16) by a Lagrange multiplier, 
A, and subtract the set from equation (15). Noting that the N weighting functions can be written 


as, 


W k (x) 


3q* (x) 
0C k (x) 


(17) 


integrating by parts the differential equation term involving the elliptic operator, and identifying 
Xa 2 with the generalized diffusional transport coefficient, k, Table 1, achieves a cancellation of 
terms. The resultant ordinary differential equation system for determination of the remaining un- 
known expansion coefficients, equation (14), is 


- /V 2 (?? )-kV 2 q*dr - plJ** V 2 q* dr 

-r 3C k j, dC k 


3 q _tt* 

k 


/■ 


3q_* 

ac. 


f* dr 


I p u* q*, dr + / 

X 3R 


ac. 


a 3 da 


/ ~ 9* da=0 

3R 3C k 


(18) 


Throughout equation (18), the superscript star notation implies an approximation to the function 
and/or parameter, consistant with equation (1 1). 

Application of these concepts on a local basis may be termed Numerical Method of Weighted 
Residuals. Establishment of an assembly procedure for generation of the global algorithm provides 
the theoretical basis for a finite element numerical solution algorithm. Since equation (1 1 ) is valid 
at a point, it is similarly valid in arbitrary subdomains, R m , of R with closures 9R m . In actuality, 
for fluid mechanics, the integral formulation transforming system behavior to a subdomain (control 
volume) can be viewed as fundamental, the more familiar differential equations having been derived 
therefrom. Hence, in each of M, disjoint interior subdomains, R m , defined by (y,z,x)eR m x[x 0 ,<»), 
and where UR m =R, define a local approximation function 

N 

Qm W = 2C mk (x)0 k (y,z) (19) 

k 


where the subscript, m, constrains the domain of equation (19) to R m . Evaluate equation (18) 
within each R m , so selected that k m (x) is adequately represented by a local value k m (x) within R m . 
This yields an ordinary differential equation system for solution of the local approximation function, 
equation (19), within R m . The global solution algorithm is established employing Boolean algebra to 
assemble the MxN equations (18) into an equation system which enforces weighted average global 
adherence to the original differential equation system, equations (1 1) through (13). 


The globally assembled equation system is similarly a first order, ordinary differential equa- 
tion system, written on the totality of local expansion coefficients, C mk (x), equation (19), that do 
not coincide with locations on the global closure, 3R, where fixed boundary conditions are applied, 
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TABLE 2 

COEFFICIENTS IN GENERAL BOUNDARY 
CONDITION STATEMENT 



Coefficients 

Boundary Condition 

3 1 

a 2 

a 3 

No-Slip at Wall 

1 

0 

0 

Slip at Wall 
Mass Injection 
Adiabatic Wall 
Specified Heat Flux 

t 

-1 

0 

0 

1 

f / 

0 

0 

1 

1 

?/ 

Temperature Dependent Flux 

t 

t 

A 

Symmetry Condition 

0 

1 

/ 0 


/ 


\ / 

t Non-zer<^ values, user specified to enforce desired 

conditions. / 

\ / 

Solution of the differential equation system, equations (4) through (7) is thus reduced to 
establishment of a numerical solution algorithm for equations/(l 1) through (13). The finite element 
method is suggested, as has been presented\by Baker (ref. 9,^10). Briefly, solution is sought to 
equation (11) using Galerkin criteria withinthe Method ofAVeighted Residuals (MWR), see Finlay- 
son and Scriven (ref. 1 1 ) by approximation of\a (each) dependent variable within the solution 
domain, R, by a series expansion of the form. 


* -► 

q (x) = 


2 C k (x) 0 k (y,z) 
k = 1 



(14) 


In equation (14), the C k (x) are unknown expansion coefficients modifying members, 0 k (y,z), of 
a set of functions that are complete in the space^of the elliptic domain. Certain of the expansion 
coefficients are evaluated by requiring that theMpproximj^e solution, equation (14), satisfies ap- 
propriate boundary conditions, equation (13)f The remaining unknown coefficients are then de- 
termined by setting to zero the differential equation residuals, formed by substitution of q* into 
equation (1 1), multiplication by a set of weighting functions, \v k , typically identical to the spatial 
approximation functions, <£ k , and integrating over the elliptic dqmain R. This produces the N 
ordinary differential equations. . / 




W R (x) L (q*) dT = Of k=l ,2,. . ,,N 


(15) 


where N corresponds to the total?number of remaining unknown coefficients in the definition of 
the approximation function, equation (14), and L(q*) is equation (11) written on the approximate 
solution. 


\ 


When non-vanishing*gradient boundary conditions exist, a generalization to MWR relaxes 
the constraints formed b^he boundary condition statement. Form the N boundary residuals, de- 
fined from equation (13/as 


/ 

3R 


W k -(I) [q*- - (a 3 -a 2 Vq*. n] da = 0 k=l,2,:..,N 



(16) 


\ 
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i.e., a 2 =0, equation (13). The order of this equation system is less than MxN by both connectivity 
of the subdomains and enforcement of fixed boundary conditions. 

The'last step to establishment of the finite element solution algorithm is selection of the 
functionals 4>^(y,z), equation (19), and specification of the finite element subdomain geometry. 
From the vast choice which exists, experience with general purpose computer codes indicates that 
many practical problems are amenable to solution using linear functionals. The natural two-dimen- 
sional finite element shape, illustrated in figure 3, is a triangle with (at least) vertex node points. 
Equation (19) constrained to a linear representation over the triangular shaped finite element do- 
main, where j } denotes a column matrix, becomes 

<4(x) = C lm (x) + C 2m (x)y + C 3m (x) z 

= {cw^jx} (20) 


v 



Figure 3. Triangular Finite Element 


Evaluating equation (20), at the three vertex nodes of the finite element, produces a linear alge- 
braic equation system which can be solved analytically for matrix elements pf|C(x)} m , Hence, an 
alternate and useful expression for q^ is 

q Mx,y,z) = I X } T [F] m {Q (x)f m (21) 

The elements of the 3x3 coefficient matrix, [n m , are known constants, derived strictly from geo- 
metrical considerations, as 





“i 

“j 

a k 

Pi 

Pj 

0k 

7j 

7j 

?k 


( 22 ) 
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where 


<*p = y q z r-y r z q 
h = z q ■ z r 


( 23 ) 


y n =y r -y , 


and the indices (p.,q,r) permute cyclicly in the order (i j,k). The elements of| x J are ( 1 ,y,z) , and the 
three elements of |Q (x)| m are the local values of the dependent variable, q^, at each node point 
location of the finite element. In equation (22), A m is the plane area of the finite element. 

The establishment of equation (21) allows direct evaluation of terms within the finite ele- 
ment algorithm, which is equation (18) constrained to the subdomain R m . The weighting functions 
for Galerkin criteria are now expressible as 


9q 


m - 


9q 


m 


ac mk 3{Q<x>! 


= in I W 

m * 


m 


(24) 


Substituting equations (21 ) and (24) into equation (18), the solution algorithm for the typical 
initial-boundary value differential equation system for the three-dimensional flow field becomes 

ICljQl m = -2 IK'lmlQlm + S Mm • (25) 

1=1 1=1 

In equation (25), the superscript prime denotes the ordinary derivative with respect to the x co- 
ordinate. The indicated finite element matrices, common to all q, are 

tci m = mj; / )x(|x| T iri m pu| m |x[ T dr ir] m 

R m 

l K "m - tn m / Jx}-k m y, { x ( T 
R m 

[K2l m -in If )x||x| T iri m (|pV| m j+|pW( m k) • 

R m 

Vj {x} T dr[r] m 
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i.e., a 2 =0, equation ( 1 3), The order of this equation system is less than MxN by both connectivity 
of the subdomains and enforcement of fixed boundary conditions. / 


The last step to establishment of the finite element solution algorithm is selection of the 
functionals equation (19), and specification of the finite element subdomain geometry. 

From th^ vast choice which exists, experience with general purpose computer codes indicates that 
many practical problems are amenable to solution using linear functionals. The^natural two-dimen- 
sional fmite.,element shape, illustrated in figure 3, is a triangle with (at least) vertex node points. 
Equation (19) constrained to a linear representation over the triangular shaped finite element do- 
main, where j ^denotes a column matrix, becomes ^ 


C lmW + C 2mWy + c 3mW z 

IcmiVi 




■ - / ; 

Figure/3. Triangular Finite Element 


Evaluating equation (20), at the^three vertex nodes of the finite element, produces a linear alge- 
braic equation system which can be solved analytically for matrix elements of |C(x)} m . Hence, an 


alternate and useful expression for q* is 


q* (x,y,z) 
m 



|x| T [ri m |Q(x)(, 


The elements of the 3x3 coefficient matrix, [P] m , are knownsconstants, derived strictly from geo- 
metrical considerations, as \\ 


/ ! a i a j a k 

ir > = 2l Pi ^ 

y M m 

7j Tj 7k 
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The source term finite element matrices, for each particular identification of q, are 

q = Y i :|L2( m = in* / 

R m 

q = u: |L2} m = -[1^ / {x}p, x dr 

R m 

q = w: |L2} m s-in ^ f |x}p !z dr 

R m 

q = H: jL2) m s ID* / jx|V(l-Pr)^()Ut* ID* H % I* 

m 

l r ! m }ul m in *W V 2 {x( T mjwl^dr 

1 l3 L -in * / |x|y 2 .(i-u)^ r)h‘|* in * jx| 

Rm 

7 2|x| T (n m {Yi| m dr (27) 

In equation (27), the pressure gradients are assumed specified functions of x and z, and the matrix 
elements of jh 1 } m are the static enthalpy at nodes of R m ,as defined by equations (8) and (9). 

Equations (25) through (27) define the finite element solution algorithm for the typical 
equation over each finite element. The MxN equations (25) are assembled, using Boolean algebra, 
into the global system which enforces a weighted average adherence to the original partial differential 
equation. The appearance of the global equation is identical to equation (25) with the subscripts m 
removed. This equation system is eligible for integration, as a system of ordinary differential equa- 
tions, written on the global column matrix, |Q(x)} , by any algorithm. The choice between implicit 
or explicit essentially reduces to whether or not elements of the [C] matrix, and the jQ(x)J vector 
on the right side of equation (25), are evaluated at the current or the next advanced x station. 

An explicit algorithm, (ref. 12), has been used to generate the reported results. For this 
choice, the elements of the j Q(x) } vector are ordered such that those nodes occurring on 3R, where 
fixed boundary conditions are applied, are loaded into the bottom partition of {Q(x)j . The sym- 
metric [Cl matrix is partitioned, correspondingly, and the rank (r) of the upper partition equals 
the number of unknowns in the { Q(x)} vector. The lower partition is contracted with the correspond- 
ing known (fixed) elements of jQ(x) } , and subtracted from both sides of equation (25). Premultiplying 
equation (25) by the inverse of the reduced [C] matrix yields, as the final global solution form for 
the finite element algorithm, the explicit equation. 

jQ)' r = [C]' 1 ^[KI] r |Q} r + 2 { LI | r ^ (28) 

The subscripts r in equation (28) indicate that the system is of the reduced rank r, and the 
contributions stemming from the lower partition operations have been lumped into an additional 
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| LI | vector. The solution of equation (28) yields ..the values of the dependent variable approxima- 
tion function q , equation (14), at the nodes of the finite element discretization of the solution do- 
main R, including those values at nodes existing on 9R where gradient boundary conditions have 
been enforced. " 

The finite element solution algorithm for global continuity follows a similar development. 
Equation (3) is strictly initial value on pv, as a function of y, with x and z appearing as parameters 
in the form of the corresponding derivatives. Hence, the finite element approximation q m , equation 
(21 ), need span only the transverse coordinate direction, as 

qj, = ji,y,y 2 ,- • -,y n l T [r] m {Q (x,z)| m (29) 

and the elements of{Q(x,z)| m are corresponding values (and derivatives) of pv at nodes of the dis- 
cretization. These nodal values are functions of both x and z; it is thus required to solve global con- 
tinuity at successive z locations, for each x station. 

The solution of equation (3) by MWR follows the conventional procedure. Identifying the 
approximation function, q m for the terms in the equation, and selecting a set of weighting func- 
tions, W k , form the weighted residual of equation (3), and integrate over the finite element domain, 
R m . This yields 

y l 

I W k (y)L(q^)dy = 0 k=l,2,...,N (30) 

y o 

which is solvable on a subdomain (finite element) basis, since equation (3) is an initial value prob- 
lem. 


Concerning the formation of particular components of equation (30), since the finite ele- 
ment approximation spans the x plane, (pw), z is known from equation (21) written on pw. How- 
ever, (pu), x is not available, since this coordinate direction is spanned by finite difference integra- 
tion. Actually, no streamwise derivatives, equations (3) through (7), are available until the flow 
field is there known. Since pv is required to evaluate these derivatives, an estimation is required 
of the (n+l) st value of (pu), x . In keeping with the use of first order accurate methods, and real- 
izing that the required variable is, in the discretized solution, actually equal to (p U)' a first-order 
accurate finite difference approximation formula is adequate, and is of the form 

(puy n+l » -2- [(pu^, .(pU) ] -(puy (o on 

x n+l' x n 

The truncation error in equation (31) equals 0.1667 h 2 (pU)'" (£), where h is the x direction step 
size, and x n < £< x n+l • 
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The source term finite element matrices, for each particular identification of q, are 


q = Y i :|L2} m = ID 


T 


/ 


m R 


{ x}S* dr 


m 


u: \ L2 \m = -t r lI I l x }P. x dr 


m 


R 


m 


q = \lL2} m =-tr ^ f {x}p, z dr 
\ R m 

<5 = H- {\n s l r <l / 

/ 

W|u| m + |w( T m in T J'\yW T lrl rJ w Lj dr 

R. 



W in m {Yi} m dr 


(27) 


In equation (27), the pressure gradients are assumed/specified functions of x and z, and the matrix 
elements of |h*| m are the static enthalpy at nodes^of R m , as defined by equations (8) and (9). 

Equations (25) through (27) define the/finite element solution algorithm for the typical 
equation over each finite element. The MxN equations (25) are assembled, using Boolean algebra, 
into the global system which enforces a weighted average adherence to the original partial differential 
equation. The appearance of the global equation is identical to equation (25) with the subscripts m 
removed. This equation system is eligibj/ for integration, as a system of ordinary differential equa- 
tions, written on the global column matrix, {Q(x) } , by any^ algorithm. The choice between implicit 
or explicit essentially reduces to whether or not elements of^the [C] matrix, and the (Q(x)} vector 
on the right side of equation (25), arc evaluated at the current or the next advanced x station. 

/ \ 

An explicit algorithm, (rbf . 12), has been used to generate the reported results. For this 
choice, the elements of the j Q(x) } vector are ordered such that those nodes occurring on 3R, where 
fixed boundary conditions aie applied, are loaded into the bottormpartition of jQ(x)} . The sym- 
metric [C] matrix is partitioned, correspondingly, and the rank (r) of the upper partition equals 
the number of unknownsfin the |Q(x)[ vector. The lower partition ^contracted with the correspond- 
ing known (fixed) elem/nts of (Q(x) } , and subtracted from both sides of equation (25). Premultiplying 
equation (25) by the^inverse of the reduced [C] matrix yields, as the final global solution form for 
the finite element algo 



Igorithm, the explicit equation. 
-1 


= [C] 


2 [KI] r j Q | 
I 


2 

I 


Li! 


\ 


( 28 ) 


\ 


The subscripts r in equation (28) indicate that the system is of the reduced rank r, and the 
contributions stemming from the lower partition operations have been lumped into, an additional 
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COMOC COMPUTER PROGRAM 


The COMOC (Computational Continuum Mechanics) general purpose finite element com- 
puter program system is coded entirely in terms of generalized non-dimensional independent and 
dependent variables. In addition to producing the results to be discussed, it has been exercised 
for problems in transient heat conduction (ref. 1 3) and the two-dimensional Navier-Stokes equa- 
tions (ref. 9, 14, and 15). 

The program consists of four basic Modules, as illustrated in figure 4. In the INPUT Mod- 
ule, the desired discretization is formed by specification of the plane coordinates (y^, z^) of the 
triad of vertex nodes for each finite element. The node numbering sequence is completely ar- 
bitrary, as is the selection of a particular discretization, and both may be chosen for output con- 
venience. In particular, nodes are first placed where output information is desired. The remain- 
der of the solution domain is then discretized to avoid finite elements with aspect ratios greater 
than about 100:1. The closure of the solution domain is recognized as the counterclockwise 
connection of the first N nodes, with N an input parameter, These node numbers are then auto- 
matically loaded into the bottom partition of the global unknown vectors, | Q (x) } . 



Figure 4. COMOC Computer Program Organization 

The appropriate non-vanishing element boundary condition coefficients a m j, equation 
(13), are included as element information for each dependent variable to be solved, since in 
general the different Q’s have individual specifications. The input phase is completed by read- 
ing non-dimensionalizing factors, various control parameters, information regarding initial con- 
ditions for the dependent variable vectors, j Q (x Q ) } , and the desired output control parameters. 

The next two Modules, figure 4, are basically DO loops on the finite elements of the dis- 
cretization. In the GEOMETRY Module, the element coefficient matrix, [T] m , is formed. 
Sequential passes through two subroutines generate and assemble the required moment distri- 
butions over areas and lines, and establish the various required finite element matrices, equations 
(26) and (27). 

The INTEGRATION Module embodies the finite difference integration algorithm for the 
system of ordinary differential equations. The basic element operation is formation of equa- 
tion (25) at a given x station, and the assembly of the element matrices into the global represent- 
ation, equation (28). In this procedure, the required order of equation (28) is automatically de- 
termined, and the matrix [C]” 1 obtained. At user-selected points, the OUTPUT Module is called 
to record the arrays within the dependent variable vectors { Q (x) ( and other desired parameters. 
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The integration algorithm presently used by COMOC is an explicit, single-step, multi-stage 
finite difference procedure with a large region of absolute stability (ref. 12). Since it is single- 
step, the changing of step-size during solution is easily accomplished. The algorithm automatically 
adjusts current step-size by comparison of a truncation error coefficient to a user-specified para- 
meter related to an allowable local error magnitude. An extensive discussion on the operation of 
this algorithm is presented by Baker and Manhardt (ref. 1 3). 


RESULTS AND DISCUSSION 


The Mechanics and Thermodynamics of Binary, Isoenergetic 
Turbulent Boundary Layer Mixing 

A considerable simplification to the general equation system, equations (3) through (7), 
can be achieved for analysis of the specific problem of interest, namely two-component, non- 
reacting cold flow rpixing of hydrogen injected into a turbulent supersonic air boundary layer 
flow on a flat plate (ref. 5 and 6). Pressure gradients in the external flow field are negligible and 
only one species continuity equation is required. As a first approximation, the flow may be assumed 
isoenegetic, and the species continuity and energy differential equation descriptions are identical. 
Hence, their respective solutions can be linearly related as 


H-H^ 

Hq-H^ 


= K 


Y-Y n 


Y -Y 
o 00 


(32) 


In equation (32), subscript zero denotes reference conditions at the point of injection, while sub- 
script infinity refers to free stream reference values. For the specific case, Y,^ vanishes and Y Q 
equals unity. The mixture specific heat is a linear function of mass fraction; for this problem, the 
specific heat is temperature independent to within about ± 5%, in the range 1 00-300° K. This variation 
was assumed negligible, and direct solution of equation (32) for static temperature of the mixture 
at any point in the solution domain is 


T s (x) 



JL u . 2 

2gJ C Pa 


1+Y (r-1) 


(33) 


Equation (33) is written on directly measurable reference values of air and hydrogen stagnation 

temperature at injection (T„ and T n ), and the computed local values of longitudinal velocity com- 

u a u h 

ponent (U) and hydrogen mass fraction (Y). The parameters g and J have their usual interpretation, 
cp a is the specific heat of air, and r is the (constant) ratio of specific heat of hydrogen to air at the air 
stagnation temperature. This assumption can be readily replaced by finite element solution of the energy 
equation (7), with temperature dependent thermophysical problems, for other problems, e.g., reacting flows 


The equation of state of the two-component mixture, equation (10), is readily reduced to 
the following expression for density, p, at any point in the solution domain assuming a constant 
static pressure distribution. 


P (x) 


Poo 



T s [1+Y (m-1)] 


(34) 
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COMOC COMPUTER PROGRAM 


The COMOC (Computational Continuum Mechanics) general purpose finite element com- 
puter program system is coded entirely in terms of generalized non-dimensional independent and 
dependent variables. In addition to producing the results to be discussed, it has been exercised 
for problems in transient heat conduction (ref. 13) and the two-dimensional Navier-Stokes equa- 
tions (ref. 9, A, and 15). / 

The program consists of four basic Modules, as illustrated in figure ^In the INPUT Mod- 
ule, the desired discretization is formed by specification of the plane coordinates (y^, z^) of the 
triad of vertex nodes ipr each finite element. The node numbering sequence is completely ar- 
bitrary, as is the selection of a particular discretization, and both may^be chosen for output con- 
venience. In particular, nodes are first placed where output information is desired. The remain- 
der of the solution domain,is then discretized to avoid finite elements with aspect ratios greater 
than about 100: 1. The closure of the solution domain is recognized as the counterclockwise 
connection of the first N nod'es, with N an input parameter. These node numbers are then auto- 
matically loaded into the bottom partition of the global unknown vectors, j Q (x)} . 


Geometry / 

\\ — LA 

1 .. ... I , / ... 

V / 

Integration Algorithm 


\ 1 / 


'Output 


Figure 4. COMOC Computer Program Organization 
/ \ 

The appropriate non-vanishing element boundary condition coefficients a m j, equation 
(13), are included as element information for each dependent variable to be solved, since in 
general the different Q’s have individua/specifications.\rhe input phase is completed by read- 

laramete 


mg non-dimensionalizing factors, various control parameters, information regarding initial con- 
ditions for the dependent variable vectors, |Q (x Q )} , and the desired output control parameters. 

./ \ 

The next two Modules, figure 4, are basically DO loops on the finite elements of the dis- 
cretization. In the GEOMETR\/Module, the element coefficient matrix, [r] m , is formed. 
Sequential passes through two^ubroutines generate and assemble the required moment distri- 
butions over areas and lines, and establish the various required finite element matrices, equations 
(26) and (27). / \ 


The INTEGRATION Module embodies the finite difference integration algorithm for the 
system of ordinary differential equations. The basic element operation's formation of equa- 
tion (25) at a given x station, and the assembly of the element matrices into the global represent- 
ation, equation (28)y^In this procedure, the required order of equation (28^ is automatically de- 
termined, and the matrix [C] r _1 obtained. At user-selected points, the OUTPUT Module is called 
to record the arrays within the dependent variable vectors { Q (x) } and other desired parameters. 
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In equation (34), T„ is the freestream (air) static temperature, T„ comes from equation (33), and 

b OO 5 

m is the ratio of air to hydrogen molecular weights. 

Finally, for the present problem of flow over a flat plate, w vanishes both at the wall and 
in the freestream, as well as along the symmetry plane bisecting the injection orifice and parallel to 
the predominant flow vector. A non-vanishing component may well exist on the remaining segment 
of the elliptic domain closure. However, since no experimental measurements are available, no 
plausible alternative exists but to assume w vanishes over the entire closure. The resultant solution 
to equation (6) is that w vanishes everywhere. 

Therefore, the partial differential equation system requiring solution for the subject cold 
flow binary mixing problem reduces to equations (3) through (5), less the respective source terms, 
in combination with equations (33) and (34). The theory for turbulent transport of mass and 
momentum is described in the next section, and the procedure for computation of molecular vis- 
cosity for the binary mixture is described in Appendix C. 

Model for Turbulent Transport Coefficient Distribution 

To date, no work is reported on the modeling of three-dimensional boundary layer flows 
for the conditions of interest in this investigation (ref. 5, 6). For example, only the incompres- 
sible case has been analytically modeled, e.g., Bradshaw (ref. 16), Nash (ref. 17, 18). Therefore, 
the most direct approach was to develop a prototype three-dimensional model by extending 
techniques proven successful for planar and axisymmetric flows. 

The three dimensional model developed herein reflects, ( 1 ) mass flux differences between 
the main flow and the jet, and (2) the turbulence due to the presence of the wall. Above the jet 
region, i.e., the region bounded by the wall and the zero concentration contour, (fig. 5) of the 
concentration profile, mixing is due to the differences in mass flux. Outside the jet and near the 
wall (large values of I z I), the turbulence is due solely to boundary layer phenomena. Within the 
region of the jet both mechanisms are active. 

Shown in figure 6 is an axisymmetric flow configuration which was studied for turbulent 
transport characterization by Morgenthaler (ref. 19), and is being further studied.* The two 
flow geometries (figs. 5 and 6) are similar in that they each consider the normal sonic injection 
of hydrogen into a supersonic air (or nitrogen) stream. Two important differences are three- 
dimensional versus axisymmetric flow, and that pressure gradients may exert influence on the 
ducted axisymmetric flow field development. Calculations have since shown that neglecting 
pressure gradient in the ducted flow case does not significantly affect the predicted hydrogen 
concentration field. Therefore, the primary difference between the flows shown in figures 5 and 6 
is the dimensionality. 

Development of the prototype three-dimensional mixing model employed the results of 
analysis for the axisymmetric configuration. For this flow detailed experimental data were avail- 
able from which distributions of eddy diffusivity coefficients were calculated using methods 
reported in reference 19. The objective was to determine whether the determined eddy diffusivity 
could be modeled using the mass defect concept and mixing length theory; i.e., techniques proven 
successful in modeling planar boundary layers (ref. 20) and axisymmetric free shear layers (ref. 

21 , 22 ). 

♦Zelazny, S.W.: Modeling of the Eddy Diffusivity of Mass in Supersonic, Axisymmetric, Ducted, 
Turbulent Flow, IOM 665:72:0320-1 :SWZ, Mar. 1972, Bell Aerospace Co., Buffalo, N. Y. 
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Figure 5. Structure of Three-Dimensional Flow Field and Mixing Region 

The mixing region is divided into an inner and outer sub-region as shown in figure 7. In 
the inner region, mixing length theory is assumed valid: hence, the eddy diffusivity is of the form 

, au 

b D . - L'-gj y I — iSc-r o < y < y d (35) 

111 0y 1 u 

where v is displacement from the wall. y d is height of inner region (defined as that value of y 

where Ep> > F.p> ), U is the mean axial velocity, and C is the mixing length, given by 
u \n 1 out 

t! - cy (36) 
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In equation (34), T c is the freestream (air) static temperature, T- comes from equation (33), and 

\ “oo 

m is theyratio of air to hydrogen molecular weights. 

Finally, for the present problem of flow over a flat plate, w vanishes both at the wall and 
in the freesrmam, as well as along the symmetry plane bisecting the injection orifice/and parallel to 
the predominant flow vector. A non-vanishing component may well exist on the remaining segment 
of the elliptic domain closure. However, since no experimental measurements are available, no 
plausible alternative exists but to assume w vanishes over the entire closure. The resultant solution 
to equation (6) is'fcjiat w vanishes everywhere. * 


Therefore, thevpartial differential equation system requiring solution for the subject cold 
flow binary mixing problem reduces to equations (3) through (5), less the respective source terms, 
in combination with equations (33) and (34). The theory for turbulent^ transport of mass and 
momentum is described irnthe next section, and the procedure for computation of molecular vis- 
cosity for the binary mixtures described in Appendix C. / 


Model for Turoulent Transport Coefficient Distribution 

\ / 

To date, no work is reported^on the modeling of three-dimensional boundary layer flows 
for the conditions of interest in this investigation (ref. 5, 6). For example, only the incompres- 
sible case has been analytically modeled e.g., Bradshaw (ref/l6), Nash (ref. 17, 18). Therefore, 
the most direct approach was to develops a prototype thre,e-dimensional model by extending 
techniques proven successful for planar and axisymmetric'flows. 


\ / 

i \ j i • c 


The three dimensional model developed hereimreflects, (1) mass flux differences between 
the main flow and the jet, and (2) the turbulence du^/to the presence of the wall. Above the jet 
region, i.e., the region bounded by the wall and tl^e # zero concentration contour, (fig. 5) of the 
concentration profile, mixing is due to the differences in mass flux. Outside the jet and near the 
wall (large values of I z I), the turbulence is due solely, to boundary layer phenomena. Within the 


region of the jet both mechanisms are active. 


L 
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Shown in figure 6 is an axisymmetric flow configuration which was studied for turbulent 
transport characterization by Morgenthaler (ref. 19), and is being further studied.* The two 
flow geometries (figs. 5 and 6) are similar irathat they eachc^nsider the normal sonic injection 
of hydrogen into a supersonic air (or nitrogen) stream. Two important differences are three- 
dimensional versus axisymmetric flow, arid that pressure gradients may exert influence on the 
ducted axisymmetric flow field development. Calculations have\ince shown that neglecting 

if ^ 

pressure gradient in the ducted flow case does not significantly affect the predicted hydrogen 


concentration field. Therefore, the primary difference between the^ 
is the dimensionality. 


lows shown in figures 5 and 6 


Development of the prototype three-dimensional mixing modeFemployed the results of 
analysis for the axisymmetric configuration. For this flow detailed experimental data were avail- 
able from which distributions of eddy diffusivity coefficients were calculated using methods 
reported in reference 19. The^objective was to determine whether the determined eddy diffusivity 
could be modeled using the mass defect concept and mixing length theory; i.e techniques proven 
successful in modeling planar boundary layers (ref. 20) and axisymmetric free'shear layers (ref. 


21 , 22 ). 


\ 


V 


*Zelazny, S.W.: Modeling of the Eddy Diffusivity of Mass in Supersonic, Axisymmetric, Ducted, 
Turbulent Flow, IOM 665:72:0320-1 :SWZ, Mar. 1972, Bell Aerospace Co., Buffalo^N. Y. 
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Figure 6. Schematic of Axisymmetric Ducted Supersonic Mixing Experimental Apparatus 
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Figure 7 Typical Velocity and Mass Fraction Profiles Downstream of Sonic Injection of H 2 
from a Circumferential Wall Slot (See Fig. 6) 


The constant c, equation (36), has been found from boundary layer data to be equal to 0.4, 
whereas the turbulent Schmidt number was assumed equal to 0.7, consistent with both boundary 
layer and free shear layer data. 

von Driest’s (ref. 23) damping factor, co, requires the eddy diffusivity to vanish at the 
wall, and is given by 

oj = (1 -e-y /A ) 2 (37) 
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where A is an empirically determined length scale expressed in terms of the wall values of shear 
stress t w and density p w . 

A = 26pA/r w /p w (38) 


The intermittency factor, 7, has been empirically modeled in a number of ways (ref. 24). 
In this study, it was found that 


1 

1 +.01f 9 


(39) 


gave a fair representation of the data where f = y/5y, and 5y is the value of y where the hydrogen 
mass fraction is one-half its wall value. 


In the outer region, the eddy diffusivity of mass is assumed directly proportional to the 
mass defect (or excess) of the concentration layer and inversely proportional to the width of the 
concentration layer; hence 

p R w 

pE D . «-/ |pU-p c U c lrdr/L (40) 

out j r^y 

where R w is the wall radius, r oy is the distance from the centerline to the point where Y equals 
0.05 Yyy^LL’ Pc * s t * ie density on the centerline, U c is the axial velocity on the centerline, Y 
is the mass fraction of hydrogen, and Y^^LL is the mass fraction of hydrogen at the wall. Also, 

L = R w -r oy (41) 


Accounting for intermittency at the outer region, the complete expression for Ep is given by 


f R w 
•' r oy 




K 7 


Ip U-p c U c I rdr 


pL 


y d <y < R 


w 


(42) 


where K is an empirical constant determined from analysis of data and found to equal 0.0072 
for a case with reliable mass balances (rhpj^ = 0.06 lbm/sec, 90° injection from a 0.03 in. 

circumferential wall slot). Thus the Ep model was completely specified in the inner and outer 
regions. The eddy diffusivity of momentum was defined by the relation, 


E M s (E D )(Sc T ) (43) 

Comparisons of the predicted Ep profiles with the experimental values led to the fol- 
lowing conclusion for the wall slot injector geometry: 

(1) The eddy diffusivities of mass, Ep, in the region away from the wall (approxi- 
mately 80 percent of the duct radius) are predicted with acceptable accuracy 
(maximum difference of 28 percent) using the mass defect as a correlating 
parameter. 
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(2) Results were inconclusive regarding suitability of mixing length theory in the 
inner region to model Ep. Additional analysis can determine whether the dif- 
ferences between predicted (using mixing length theory) and experimental 
Ep’s have a significant effect on the predicted hydrogen concentration. 

For calculating the transport coefficients in the outer region of the three-dimensional flow 
field, equation (42) was replaced by equation (44). 

E d = KYl/pL' (44) 

where K' is an empirical constant determined from data analysis, y' is the intermitancy factor, equation 
(39), but with an argument of y/5, L' is a characteristic length of the mixing region defined as the half 
height of the mixing layer on the centerplane, 6, and I is the mass defect in the three-dimensional field. 
The evaluation of mass defect in the three-dimensional flow field requires integration over the elliptic 
solution domain of the form 

I(x) =f R IpV-p^U^ldr (45) 

where the subscripts infinity refer to local free stream reference values. Within the concepts of 
finite element formalisms, an entirely equivalent form, using the appropriate approximation func- 
tion description, equation (19), and summation over the finite elements of the discretization is 


I (x) =ZI m (x) 

m 111 



IP Um - P„ U„ |dydz 


(46) 


The form of pU^is known and PooU^ is a local scalar constant. Hence, using equation (21), the 
mass defect computation becomes 

T f (xi 

Jd ) X ( 


I (x) = 2 l m (x) = 2 


m 


m 


[rhou! t in 


)m 


m J R m V 


dr - Poo U«,A m 


(47) 


where A m is the plane area of the m tla finite element. Since both |rHOu] m and p x may be 
functions of x, evaluation of equation (47) is required at each longitudinal computational station. 
Shown in Figure 8 is a typical COMOC computed turbulent mass mixing coefficient, Ep, distribution 
through the boundary layer, using equations (35) and (44) for the multiple-jet configuration of Rogers 
(ref. 6). | 


The turbulent mixing coefficient distribution used for prediction of momentum transport, E m , 
at any point in the flow field requires specification of a turbulent Schmidt number, Scj. As with all 
transport coefficients, COMOC can accept multiple methods of specification for such distributions. 
Initial studies employed a uniform constant Scj, equal to 0.7, which neglected the three-dimensionality 
of the flow field and consistently predicted over-diffusion of momentum near the wall. Considerable 
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Figure 8. Computed Turbulent Mass Transfer Coefficient Profile on Centerplane (z/D = 0.0) of 

Multiple-Jet Geometry, q r = 1 .0, x/D = 1 20 


experimentation led to an empirical relationship for a geometric Schmidt number distribution, that 
enhanced agreement with data, of the form 


r(0 ) 

Sc-p = 1. - R(0) r < R (48) 

In equation (48). R(0) describes an ellipse with major axis parallel to the plate and normal to the' 
predominant flow direction, whose center is located at the approximate centroid of the imbedded 
hydrogen jet (fig. 2). The ellipse selected for the present studies was specified independent of 
downstream station and established by equation (49). 


y/D - 2.5 
2.5 



(49) 
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Results were inconclusive regarding suitability of mixing length theory in/the 
inner region to model Ep. Additional analysis can determine whether the dif- 
ferences between predicted (using mixing length theory) and experimental 
Ep’s have a significant effect on the predicted hydrogen concentration. 

/ 

For calculating the transport coefficients in the outer region of the three-dimensional flow 
field, equation\>42) was replaced by equation (44). 



E D = 



'y'l/pV 



(44) 


where K' is an empirical constant determined from data analysis, 7 ' is the intermitancy factor, equation 
(39), but with an argument of y/5, L' is a characteristic length of the mixing region defined as the half 
height of the mixing layer\on the centerplane, 5, and I is the mass defect in the three-dimensional field. 
The evaluation of mass defect in the three-dimensional flow field requires integration over the elliptic 
solution domain of the fomv 
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where the subscripts infinity refer to local free stream reference values. Within the concepts of 
finite element formalisms, an entirely, equivalent form, using the^appropriate approximation func- 
tion description, equation (19), and summation over the finite elements of the discretization is 


I(x) =Sl m (x) = 2 


‘m v-, — / 

m m J R 


\ / 
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The form of pU^is known and p^U^ is a local scalar constant. Hence, using equation (21), the 


mass defect computation becomes 
I (x) = 2 I m (x) - Z 


m 


m 



[RHOUr IVh 
! lm V m J R r 


jxj dr - Poo Uoo A m 


(47) 


where A m is the plane area of the m ** 1 finite element. Since both I RHOul and p„ U„ may be 
functions ot x, evaluation of equation (47) is required at each longitudinal computational station. 
Shown in Figure 8 is a typical COMOC computed turbulent mass mixing coefficient, Eq, distribution 
through the boundary layer, using equations (3^5) and (44) for'fhe multiple-jet configuration of Rogers 
(ref. 6 ). / 


The turbulent mixing coefficient distribution used for prediction of momentum transport, E m , 
at any point in the flow field requires specification of a turbulent Schmidt number, Sc-p. As with all 
transport coefficients, COMOC can accept multiple methods of specification for such distributions. 
Initial studies employed a uniform constant Sc-p, equal to 0.7, which neglected the three-dimensionality 
of the flow field and consistently predicted over-diffusion of momentum^near the wall. Considerable 
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Additionally, r(0) is the magnitude of a position vector to any point within the ellipse and aligned 
parallel to R(0), and determined as 

r(0) = '/(y/D - 2.5 ) 2 + (z/D) 2 (50) 

Along the plate surface, the Schmidt number was uniformly set equal to 0.7. Elsewhere outside 
the ellipse, the Schmidt number distribution was held constant at those values used in the initial 
studies that determined equation (48). 

Prediction of Mass and Momentum Transport Following Transverse 
Injection of Hydrogen from a Single Orifice 

At some location downstream from transverse injection, as a function of the dynamic pressure 
ratio q r , boundary layer reattachment occurs, (fig. 5a), and the three-dimensional elliptic boundary 
layer equations become valid. For the configuration of Rogers (ref. 5), reattachment should occur by 
station x/D=30, where D is the injection orifice diameter. Therefore, evaluation of the prototype 
three-dimensional turbulent mixing model was performed in the downstream region 30 < x/D < 120. 

It is necessary to establish initial conditions for the planar distribution of all dependent vari- 
ables at x/D=30. The original raw data of Rogers (ref. 5) consisted of pitot-static and hydrogen sam- 
pling data points on four traverses in the (y, z) plane (fig. 2). The vertical traverse, in the streamwise 
plane through the injection orifice and parallel to the y axis, determined the y locations where hydro- 
gen concentration was maximum, (y/D) max , and where it vanished, i.e., the freestream boundary. 

Three horizontal data traverses were then made parallel to the plate surface, at y locations correspond- 
ing to (y/D) max , and at a location each side of (y/D) rn | ax approximately bisecting the respective do- 
mains. Shown in figure 9 are the horizontal traverse data at x/D=30 for q r =1.0. The curves appear of 
Gaussian shape; an apparent symmetry plane exists which is displaced from the geometric plane of 
symmetry, i.e., z=0, as Rogers observed. 

Although the entire flow field could be computed numerically, the strong appearance of a data 
symmetry plane suggests establishing a corresponding solution domain. Therefore, these data were in- 
put to a cubic spline interpolation computer program that established the z location of the data sym- 
metry plane via a minimization criteria. For the data of figure 9, the symmetry plane was located at 
z/D=-0.7, in agreement with folding the data manually and matching the “wings” of the Gaussian 
shaped curves. Having thus established the symmetry plane, the spline package curve fit both the Y 
and u displaced data, and evaluated the interpolation polynomials at node points of the finite element 
discretization of the elliptic solution domain (fig. 10). A number of different discretizations were eval- 
uated; in each case, the spline package provided appropriate initial data. 

Shown in figure 1 1 , for planes z=constant of the discretization of figure 1 0, are the spline-computed 
data fits compared to the span of the measured raw data, and the subsequent data range for the “best 
symmetry plane” determination. Data symmetry occurs to within about ±0.2% of the maximum hy- 
drogen deviation, and computations were performed using only the half-plane discretization of figure 
10. These same manipulation operations were performed on Rogers data at stations x/D=60 and 120 
to establish comparison bases. Shown in figure 12 is the corresponding, and more familiar (see ref. 5) 
contour distribution of hydrogen mass fraction at x/D=30, as computed by the symmetrized spline func- 
tion interpolation of the data. In both figures 1 1 and 12, note that the contours do not penetrate to 
the plate surface, since no data are available below the lowest horizontal data survey except on the 
geometry symmetry plane. 


21 



Station x/D = 30, From Rogers (ref. 5) 



Figure 9. Experimentally Measured Mass Fraction for Injection from a Single Orifice 


Shown in figure 13 is the COMOC computed hydrogen mass fraction contour distribution at 
station x/D=60. Superimposed for comparison purposes are the (spread of) data as spline-interpolated 
for the Pest symmetry plane fit. Agreement of the results along the symmetry plane, z = 0., is excellent. 
These results were achieved for the turbulent mixing model previously described for K ’= 0.1 (equation 44). 
Transition from mixing length to mass defect modeling occurred between 0.6 and 1 .0 injection diameters 
above the plate, across the entire pattern. The lateral spreading of the jet (parallel to z axis) is predicted 
accurately near the wall, but excessive concentration levels are predicted in the middle region of the 
pattern. 

The transition from the initial distribution, and significant detail on solution accuracy, are pre- 
sented in figure 14. which contains data and computed concentration profiles along planes z=constant 
at x/D=60. The indicated disagreement between data and computations, figure 13. may be directly 
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Additionally, r(0) is the magnitude of a position vector to any point within the ellipse and aligned 
parallel to K{0), and determined as 


r(0) = v(y/D - 2.5) 2 + (z/D) 2 


Along the plate surface, the Schmidt number was uniformly set equal to 0.7. Elsewhere outside 
the ellipse, the Schmidt number distribution was held constant at those values used in the initial 


studies that determine^ equation (48). 


Prediction of Mass and Momentum Transport Following Transverse 
\ Injection of Hydrogen from a Single Orifice/ 


At some location downstream from transverse injection, as a function of the dynamic pressure 
ratio q r , boundary layer reattachment occurs, (fig. 5a), and the threendimensional elliptic boundary 
layer equations become valid. For the configuration of Rogers (ref/5), reattachment should occur by 
station x/D=30, where D is the injection orifice diameter. Therefore, evaluation of the prototype 
three-dimensional turbulent mixing model was performed in the/iownstream region 30 < x/D < 120. 


It is necessary to establish initialconditions for the planar distribution of all dependent vari- 
ables at x/D=30. The original raw data ofyRogers (ref. 5) consisted of pitot-static and hydrogen sam- 
pling data points on four traverses in the (y. z) plane (fig. 2). The vertical traverse, in the streamwise 
plane through the injection orifice and parallel to the y axis, determined the y locations where hydro- 
gen concentration was maximum, (y/D) max , and where’it vanished, i.e., the freestream boundary. 

Three horizontal data traverses were then madeS^arallel to the plate surface, at y locations correspond- 
ing to (y/D) max , and at a location each side of (wD,)^ ax approximately bisecting the respective do- 
mains. Shown in figure 9 are the horizontal traverse data at x/D=30 for q r =l .0. The curves appear of 
Gaussian shape; an apparent symmetry plane exists which is displaced from the geometric plane of 
symmetry, i.e., z=0, as Rogers observed. / \ 

/ \ 

Although the entire flow field could be computee^numerically, the strong appearance of a data 
symmetry plane suggests establishing a corresponding solution domain. Therefore, these data were in- 
put to a cubic spline interpolation computer program that established the z location of the data sym- 
metry plane via a minimization criteria. For the data of figure 9, the symmetry plane was located at 
z/D=-0.7, in agreement with folding the/data manually and matching the “wings” of the Gaussian 
shaped curves. Having thus established' the symmetry plane, theyspline package curve fit both the Y 
and u displaced data, and evaluated the interpolation polynomial^ at node points of the finite element 
discretization of the elliptic solution/domain (fig. 10). A number o^f different discretizations were eval- 
uated; in each case, the spline package provided appropriate initial data. 

/ \ 

Shown in figure 1 1 , for planes z=constant of the discretization of figure 1 0, are the spline-computed 


data fits compared to the span^f the measured raw data, and the subsequent data range for the “best 
symmetry plane” determination. Data symmetry occurs to within aboufi±0.2% of the maximum hy- 
drogen deviation, and computations were performed using only the half-plane discretization of figure 
10. These same manipulafion operations were performed on Rogers data attestations x/D=60 and 120 
to establish comparison bases. Shown in figure 12 is the corresponding, and more familiar (see ref. 5) 
contour distribution of hydrogen mass fraction at x/D=30, as computed by the^symmetrized spline func- 
tion interpolation of the data. In both figures 1 1 and 12, note that the contour.;? do not penetrate to 
the plate surface, since no data are available below the lowest horizontal data survey except on the 
geometry symmetry- plane. \ 
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Figure 1 0. Finite Element Discretion of Elliptic Solution Domain 
for Single-Jet Injection Analyses 
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Figure 1 1 . Single-Jet Hydrogen Mass Fraction Distributions 
i at Station x/D = 30, q r = 1 .0, After Rogers (Ref. 5) 

assessed as lack of hydrogen diffusion from the core region in the lateral direction. Since the turbulent 
transport model is an extension of two-dimensional concepts, it is insensitive to variable gradients par- 
allel to the z-axis. However, the experimental data show that in this region, u,y and u, z are of about 
' equal magnitude. This is readily observed in figure 15, which is a three-dimensional surface plot of the 
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Symbols are Best Symmetry Plane Fit for Data of Rogers (Ref. 5) 
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Figure 12. Cubic Spline Interpolated Hydrogen Mass Fraction Contours for 

Single-Jet, q r = 1 .0, x/D = 30. 


Symbols are Best Symmetry Plane Fit for Data of Rogers (Ref. 5) 
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Figure 13. Comparison Between COMOC Computed Hydrogen Mass Fraction Contours 
and Experimental Data for Single - Jet, q r = 1 .0, x/D = 60. 

velocity distribution at x/D=30, as observed from a viewpoint beneath the plate surface. The superim- 
posed grid coincides with the finite element discretization (fig. 10) and the hydrogen jet is imbedded 
within the centroidal indentation. Obviously, these three-dimensional effects are important, and 
should form an integral part of future development studies. The lack of data and subsequent neglect 
of lateral velocity (w) would also exert an influence on spreading. 
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Figure M. Single-Jet Hydrogen Mass Fraction Distributions 
at Station x/D = 30, q r = 1 .0, After Rogers (R'ef. 5) 
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assessed as lack of hydrogemdiffusion from the core region in the lateral direction. Since the turbulent 
transport model is an extension of two-dimensional concepts, it is insensitiveUo variable gradients par- 
allel to the z-axis. However, the experimental data show that in this region, u^, and u, z are of about 
equal magnitude. Thisis readily observed in figure 15, which is a three-dimensi'onal surface plot of the 
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Symbol Description 

— Initial Conditions, x/D = 30 

A Raw Data Measurements 

▲ Best Symmetry Plane Determination 

— COMOC Computation 



Hydrogen Mass Fraction - Y (%) 

Figure 1 4. Computed Single-Jet Hydrogen Mass Fraction Distribution 
at Station x/D = 60, q r = 1.0 


Note also in figure 14, that the initial data profiles (dashed lines) have been extrapolated to the 
plate surface, as is required to establish initial conditions for hydrogen mass fraction, Y. Only after 
considerable scrutinizing of the data, and extensive numerical experimentation with boundary condi- 
tions for Y on the plate surface, were the illustrated results achieved. Referring to figure 1 1, con- 
siderable option exists for extrapolating the data at planes z=0. and z=l .0. It appears that the hydro- 
gen jet is held off the wall in this region by some mechanism, while the Y contours for larger z quite 
obviously will intersect the wall in an approximately normal fashion. Employing a vanishing gradient 
boundary condition for Y across the entire plate surface always resulted in the maximum local hydro- 
gen concentration sinking to the wall by, or shortly after, station 60. This result suggested that alter- 
nate wall boundary conditions be considered. Looking at experimental data for stations 30, 60, and 
120 (fig. 16), one might conclude that the wall concentration of hydrogen never exceeds about 2% 
directly underneath the injection pattern, and that a local maximum may well exist off-axis, i.e., z^O. 
The computational solutions, leading to the results shown in figures 14 and 16, numerically “froze” 
the concentration of hydrogen at 2% at the wall nodes, located at z/D=0. and z/D=1.0, while enforc- 
ing a vanishing normal gradient at the remaining wall nodes, i.e., z/D>1.0. For these conditions, the 
indicated good agreement near the centerplane was achieved. 

It is instructive to note that these observations, and the versatility of the numerical experi- 
mentation capability to routinely change boundary conditions and test theories, may exert consider- 
able impact upon design. In figure 5, reattachment of the boundary layer is estimated to occur at 
some point downstream of injection. This has been largely confirmed; however, the physics of the 
reattachment appear to isolate the peak hydrogen concentration from the plate. The impact of this 
phenomenon on actual combustion, and wall cooling requirements, could be rather important for actual 
combustor design. 
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Figure 15. Isometric View of Tliree-Dimensional Velocity Surface for Single-Jet Configuration, 

x/D = 30, After Rogers (Ref. 5) 

Figures 1 7 and 18 compare the computed solutions for longitudinal velocity distribution to 
data at stations x/D=60 and 1 20 as well as the initial distributions. Agreement is generally good. 

The experimental data are quite sparse; the apparent acceleration of How near the plat surface be- 
tween x/D=30 and 60 was computationally encouraged by the Schmidt number distribution described 
by equation (48). The CPU execution time on the IBM 360/65 digital computer, including all 
turbulence computations and detailed output was 1070 seconds with all operations performed using 
single precision arithmetic. 


Prediction of Mass and Momentum Transport Following 
Transverse Injection From Multiple Orifices 

The computational study described for the single-jet configuration was repeated for data of 
Rogers (ref. 6) obtained for transverse injection (at q r =l .0 into a Mach 4 airstream) from a row of 
orifices, aligned perpendicular to the main flow vector with a uniform separation distance of 12.5 
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A Raw Data Measurements 

A\ Best Symmetry Plane Determination 

- COMOC Computation 
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Figure 14. Computed Single-Jet Hydrogen Mass Fraction Distribution 
at Statfo^n x/D = 60, / r = 1 .0 

Note also in figure 14, that the initial dVa profile/ (dashed lines) have been extrapolated to the 
plate surface, as is required to establish initial coWiitionsfor hydrogen mass fraction, Y. Only after 
considerable scrutinizing of the data, and extensivk numerical experimentation with boundary condi- 
tions for Y on the plate surface, were the illustratedjresults achieved. Referring to figure 1 1, con- 
ition exists for extrapolating the data atiplanes z=0. and z=l .0. It appi 


siderable option 




appears that the hydro- 


gen jet is held off the wall in this region by some mechanism, while the Y contours for larger z quite 
obviously will intersect the wall in an approximately normal fashion. Employing a vanishing gradient 
boundary condition for Y across the entire plate surface always resulted in the maximum local hydro- 
gen concentration sinking to the wall by, or shortly after, station 60. This result suggested that alter- 
nate wall boundary conditions be considered/ Looking at experimental data for stations 30, 60, and 
120 (fig. 1 6), one might conclude that the^all concentration^ hydrogen never exceeds about 2% 
directly underneath the injection pattern, ^and that a local maximum may well exist off-axis, i.e., z#0. 
The computational solutions, leading to the results shown in figures 14 and 16, numerically “froze” 


the concentration of hydrogen at 2% atfthe wall nodes, located at'^/D=0. and z/D=1.0, while enforc- 
ing a vanishing normal gradient at the^remaining wall nodes, i.e., z/E(> 1 .0. For these conditions, the 
indicated good agreement near the centerplane was achieved. 


It is instructive to note that these observations, and the versatility of the numerical experi- 
mentation capability to routinely/ehange boundary conditions and test theories, may exert consider- 
able impact upon design. In figure 5, reattachment of the boundary layers estimated to occur at 
some point downstream of injection. This has been largely confirmed; however, the physics of the 
reattachment appear to isolate the peak hydrogen concentration from the plate. The impact of this 
phenomenon on actual combustion, and wall cooling requirements, could be father important for actual 
combustor design. 
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Figure ] 7. Computed Single-Jet Longitudinal Velocity Distribution 
at Station x/D = 60, q r = 1 
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Figure 1 8. Computed Single-Jet Longitudinal Velocity Distribution at Station x/D = 1 20, q r - 1 0 

orifice diameters (D). Discretization of the finite element domain parallel to the y axis, figure 10, 
remained the same. The discretization in the lateral direction was refined such that the node col- 
umns spanned -6.25 < z/D < 6.25, i.e., the half-width of orifice separation. 


The raw data were again analyzed using the cubic-spline interpolation package. A data sym- 
metry plane was established by the minimization procedure. Shown in figure 19 are the spline- 
interpolated data profiles on planes z=constant at station x/D=30, as well as the best symmetry plane 
data spread. Therefore, as with the single-jet study, a half-plane computational domain was suffic- 
ient. Vanishing gradient boundary conditions were then applied along both lateral planes to simulate 
multiple injectors. 

Shown in figures 20 and 2 1 are the COMOC computed hydrogen mass fraction profile dis- 
tributions at stations x/D=60 and x/D=120, and figure 22 displays the more familiar contour plot 
of the hydrogen distribution at x/D=120. These results were obtained using the identical mixing 
model of the single jet, i.e., K=0.1 (equation 44) with transition from mixing length to mass defect 
occurring in the region 0.6 < y/D <1.1. Figure 20 indicates that the center plane diffusion is some- 
what over-predicted, while a considerable improvement between data and computations has occurred 
in the lateral region; i.e., z/D > 2.0, using the same “frozen” hydrogen boundary condition at the 
plate surface underneath the jet centroid (only). These same conclusions hold at x/'D= 1 20 (fig. 21) 
and figure 22 illustrates how the contour patterns merge between jets for the multiple injector con- 
figuration. Agreement between computed and measured velocity distributions at both downstream, 
stations was good. 
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Figure 16. Computed Single-Jet Hydrogen Mass Fraction 
Distribution at Station x/D = L20, q r = 1 
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Figure/] 7. Computed Single-Jet Longitudinal Velocity Distribution 
/ at Station x/D = 60, q r = 1 \ 
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Symbol Description' 

Initial Conditions, x/D = 30, 

O Raw Data Measurements 

• Best Symmetry Plane Determination 

COMOC Computation 



Hydrogen Mass Fraction - Y (%) 


Figure 21. Computed Multi-Jet Hydrogen Mass Fraction Comparison 
Distribution at Station x/D = 120, q r = 1 


The apparent over-diffusion within the core of the multi-jet configuration may be a reflection 
of the overall lower hydrogen levels used as initial data. Shown in figure 23 are comparisons of pre- 
dicted and measured decay of the maximum local hydrogen concentration as a function of x/D for the 
single- and multiple-jet geometries. The computational curves are essentially parallel: the multiple-jet 
data are coincident with the single-jet except at x/D=30 where initial conditions were established. The 
decay of the single jet was accurately predicted by COMOC. The indicated lack of agreement for the 
multiple jet is perhaps a reflection of differences in initial conditions. 


Computational solutions can produce incredible volumes of output data, and integral parameters 
are often times useful to establish trends. Since the performance prediction of combustion systems is 
one ultimate goal of design studies. COMOC was adapted to compute the integral “mixing efficiency” 
parameter, r\ of reference 6, defined as the percentage of hydrogen that could completely react 
in a given concentration and velocity distribution within an air boundary layer. For the three- 
dimensional flow fields considered, r? is a point function of longitudinal displacement, x/D. In plane 
regions that are locally fuel-lean, all the hydrogen could potentially react. In fuel-rich regions, stoichio- 
metry limits the reaction. The “mixing efficiency” is, therefore, the sum of these terms, divided by the 
total hydrogen available, i.e.. 


f puYdydz + f puR-d-Yjdydz 

R, R., 


T? ( X ) = •- 


/ puYdydz 


R|F R , 


( 51 ) 
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Figure 20. Computed Multi-Jet Hydrogen Mass FractiorijComparison 
/ Distribution at Station x/D = 60, q r = 1 \ 
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Figure 22. Comparison between COMOC-Computed Mass Fraction Contours and Experimental Data for 

Multiple-Jet, x/D = 120, q r = 1.0 
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Figure 23. Downstream Decay of Maximum Hydrogen Concentration for 
Single and Multiple-Jet Transverse Injection 

Regions l and 11 are respectively luel-iean and tuel-rich, and R s (=0.029) is the stoichiometric mixture 
ratio for the hydrogen/air system. Note that the denominator of equation (51) represents total hydro- 
gen flow rate, evaluation of which is discussed in detail in Appendix A. Equation (51) is readily eval- 
uated as a summation over the Finite elements of the discretization, as 

M 

i?(x) 55 I V m (x) (52) 

m= I 

Determination of region 1 or II dependence is made on an averaged Y concentration basis within each 
finite element. 

COMOC evaluates equation (52) at each output station. Shown in figure 24 are computed mix- 
ing efficiency distributions for the single- and multiple-jet configuration. The multiple-jet efficiency 
lies uniformly higher than the single jet, as a direct consequence at least of overall lower hydrogen con- 
centration levels starting with the data for initial conditions (see fig. 23). The indicated 15% relative 
mixing improvement of the multiple jet by station x/D-85 may also be influenced by discretization 
and initial condition differences. The “virtual source” studies, discussed in the next section, circum- 
vent these difficulties by providing a uniform basis for relative comparison. 

Since equation (51) is a relative measure, it should be rather insensitive to computed conserva- 
tion phenomena. For both cases, the numerical evaluation of the denominator of equation (5 1 ) indi- 
cated about an 8% computational hydrogen flow loss by x/D=60, with an additional 1 0% loss accruing 
by x/D=120. As discussed in Appendix A, these apparent losses appear acceptable in terms of comput- 
ed solution accuracy for the present study. They may likely become unacceptable for assessing more 
accurate mixing model hypotheses. A finer solution domain discretization, proven to increase solution 
accuracy (see Appendix A) should then be used. 
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Figure 22. Comparison between COMOC-Computed Mass Fraction Contours and Experimental Data for 

Multiple-Jet, x/D = 1 20, q r = 1.0 



Figure 24. Computed Mixing Efficiency for Single and Multiple Transverse Injection 


Evaluation of a “Virtual Source” Concept 


The discussed results for single- and multiple-jet injection geometries indicate that the proto- 
type turbulence model has captured at least the essential character of binary mixing within this three- 
dimensional boundary layer flow. However, it is equally evident that initial conditions do exert a 
strong influence on computed results for each distinct configuration. Starting computational solutions 
from a “similarity” condition is an alternative to use of initial data, and one candidate method for the 
present study is a “virtual source”. Within this concept, the complex transverse injection phenomena 
is computationally replaced by a hydrogen jet imbedded within the air boundary layer flow, the dis- 
tinguishing features of which are solely dependent upon freestream and injectant parameters. Explora- 
tory results for a subsonic configuration, similar to the single-jet geometry, are discussed in reference 29. 


A theoretical model for establishing initial conditions for a virtual source was derived. The 
model captures the essence of the barrel shock-Mach disk hypothesis for injectant-freestream equilibra- 
tion, (e.g., ref. 25), and the details are discussed in Appendix B. For the present study, the virtual 
source was established in the plane of injection, i.e., x/D = 0.0, and was assumed to be of elliptical 
cross-sectional shape and composed of 100% hydrogen. The virtual source mass flow through the 
ellipse computationally coincided with the hydrogen mass flow rates of Rogers (ref. 5 and 6), for each 
dynamic pressure ratio q r . The velocity of the virtual source injectant was uniformly assumed equal to 
that velocity (1500 ft/sec or 457 m/sec) computed to occur immediately downstream of the Mach disk, 
and was independent of dynamic pressure ratio. Everywhere exterior to the ellipse, the hydrogen con- 
centration identically vanished. Except for directly above the ellipse, the virtual source was assumed 
imbedded within the turbulent boundary layer developing in the absence of injection. Above the ellipse, 
to account for displacement effects of the barrel shock, the air velocity was assumed initially uniform 
at the freestream value. 
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Figure 25 is a three-dimensional representation of the initial station velocity surface; the unity 
hydrogen concentration is imbedded within the centroidal depression. As before, the superimposed 
grid is the finite element mesh, only one-half of which was computationally required (from symmetry). 
Shown in figure 26 is the corresponding planar view of the discretization for the multiple-jet configura- 
tion (ref. 6). Three dynamic pressure ratios (q r = 0.5, 1.0, 1 .5) were analyzed; the ellipse corresponding 
to q r = 1.5 is shown in figure 26. For q r = 1.0, the reduced size of the injectant region excluded the en- 
circled nodes (fig. 26) while both the encircled and crossed node points were exterior to the virtual 
source ellipse for the case q r = 0.5. 



Figure 25. Initial Condition Velocity Surface for Virtual Source 
Simulation of Transverse Injection 
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Figure 24. Computed MixingvEfficiency for Single and Multiple Transverse Injection 

\ / 

Evaluation of a “Virtual Source” Concept 
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The discussed results for single- and multiple-jet injection geometries indicate that the proto- 
type turbulence model has captured at least the/essential character of binary mixing within this three- 
dimensional boundary layer flow. However, itns equally evident that initial conditions do exert a 
strong influence on computed results for each distinct configuration. Starting computational solutions 
from a “similarity” condition is an alternative to use of initial data, and one candidate method for the 
present study is a “virtual source”. Within this conceptMhe complex transverse injection phenomena 
is computationally replaced by a hydrogen jet imbedded wdthin the air boundary layer flow, the dis- 
tinguishing features of which are solely/dependent upon freestream and injectant parameters. Explora- 
tory results for a subsonic configuration, similar to the singleyet geometry, are discussed in reference 29. 


A theoretical model for establishing initial conditions fora virtual source was derived. The 
model captures the essence of the barrel shock-Mach disk hypothesis for injectant-freestream equilibra- 
tion, (e.g., ref. 25), and the details are discussed in Appendix B. FoVthe present study, the virtual 
source was established in theqdane of injection, i.e., x/D = 0.0, and Vas assumed to be of elliptical 
cross-sectional shape and cjamposed of 100% hydrogen. The virtual source mass flow through the 
ellipse computationally coincided with the hydrogen mass flow rates of Rogers (ref. 5 and 6), for each 
dynamic pressure ratio qf. The velocity of the virtual source injectant wa\ uniformly assumed equal to 
that velocity (1500 ft/sec or 457 m/sec) computed to occur immediately downstream of the Mach disk, 
and was independent^ dynamic pressure ratio. Everywhere exterior tc the\llipse, the hydrogen con- 
centration identically vanished. Except for directly above the ellipse, the virtual source was assumed 
imbedded within/he turbulent boundary layer developing in the absence of inj$gtion. Above the ellipse, 
to account for displacement effects of the barrel shock, the air velocity was assui 
at the freestream value. 
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Figure 26 Solution Domain Finite Element Discretization for Virtual Source Studies 

The turbulence transport model used for the virtual source study was identical to that previously 
described, except for minor variations. The mixing length hypothesis was uniformly enforced in the 
initial region until the minimum velocity in the virtual source depression accelerated to within ~ 2% of 
the corresponding velocity without injection. This occurred within 3 to 4 diameters downstream of 
the injection plane for all q r Downstream of x/D = 4, transition from mixing length to mass defect 
occurred between 0.5 and 1.0 diameters above the plate surface. Due to the relatively small density 
within the virtual source, the computed mass defect, equation (45), was quite large in the initial down- 
stream region. From extensive experimentation, a smaller constant (K = 0.01 , equation (44)) was found 
to be effective for all three dynamic pressure ratios studied. The turbulent Schmidt number distribution 

coincided with that of the single- and multiple-jet studies discussed. 

1 1 


Shown in figures 27 through 29 are the results of the virtual source computational simulation of 
the multiple injection configuration of ref. 6 for dynamic pressure ratios (q r ) of 0.5, 1 .0, and 1 .5, re- 
spectively. Superimposed are appropriate single-jet data from ref. 5, as well as the multiple-jet data. 

The computational results are compared, in the region of 0 < x/D < 60, on the, multiple bases of: 


(a) maximum concentration level of hydrogen at x/D 

(b) trajectory of the maximum hydrogen concentration level (i.e., transverse displacement 
above the plate of the local maximum hydrogen level) 

(c) lateral spreading of the jet as determined by the z/D coordinate of the hydrogen concentra- 
tion equal to 10% of the local maximum level. 


35 



The computed mixing efficiency parameter r?, equation (48), completes each data set. Listed in table 3 
are the parameters pertinent to each solution, the downstream integration station at which COMOC 
automatically switched from the mixing length to mass defect model within the jet region, and the total 
CPU execution time for the IBM 360/65 computer. 

TABLE 3 

PARAMETERS FOR THE VIRTUAL SOURCE COMPUTATIONS, 

M^ = 4.57 


Dynamic 
Pressure 
Ratio, q r 

Virtual Source Flow Area 
(Square Orifice Diameters) 

Downstream 
Station, x/D, 
for Core Region 
Mass Defect 3 

CPU Time 
Seconds 4 

Ref. 5 1 

COMOC 2 

0.5 

2.28 

2.22 

4.13 

916. 

1.0 

2.72 

2.94 

3.53 

943. 

1.5 

4.08 

3.48 

4.04 

1111. 


^ Initial elliptical cross sectional area required to match test injection conditions of 
reference 5. 

Numerical approximation to conditions of reference 5. 

2 Automatically established at first integration station past a user-specified value of 
x/D. 

4 Single precision on IBM 360/65 computer. 


The agreement of the computed solutions with data, for q r = 0.5 and 1 .0, is quite good (maxi- 
mum deviation ~ 20%) on all comparison bases. Poorer (40% deviation) agreement with the data 
occurs at q r = 1.5, (fig. 29) on both trajectory and lateral spreading. This dynamic pressure ratio pro- 
duces a “stiffer” jet, the influence of which is especially noticeable within the region 0 < x/D < 30, 
for the illustrated solution differences on both trajectory and maximum concentration. The lateral 
spreading disagreement may result primarily from the relatively undersized initial condition ellipse, 
(table 3). This could have been readily corrected by initial condition and discretization changes, but 
the uniform comparison basis of the computations for the three values of q r would then have been 
disturbed. 

Agreement between the COMOC virtual source computations and experimental data indi- 
cates a potential value for practical design studies. In providing such support, it is significantly im- 
portant that arbitraty parameters were set uniformly constant for each q r case, as distinct from fine 
tuning each computation to match specific data sets. It is thus possible to directly compare the 
results of computation without significant interpretation. As opposed to scrutiny of comparative 
flow field details, the valuable design criteria may well be the mixing efficiency comparison previous- 
ly shown in figure 1. Since combustion systems are of ultimate interest, and if the reaction is fast 
(hence diffusion controlled), the computation of mixing efficiency on a cold flow basis may provide 
some degree of design guidance for hardware development. These reported results are representative 
of an initial theoretical study; certainly, with the refinements that are possible in modeling and the 
significant versatility of the developing COMOC system, analytical support to design optimization 
appears practical. 
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Figure 26 Solution Domain Finite Element Discretization for Virtual Source Studies 
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The turbulence transport model used for the virtual source study was identical to that previously 
described, except for minor variations. ^Fhe mixing^length hypothesis was uniformly enforced in the 
initial region until the minimum velocity in the virtual source depression accelerated to within ~ 2% of 
the corresponding velocity without injection. This occurred within 3 to 4 diameters downstream of 
the injection plane for all q r Downstream of x/D = 4, transition from mixing length to mass defect 
occurred between 0.5 and 1.0 diameters above the plate surface. Due to the relatively small density 
within the virtual source, the computed mass defect, equation (45), was quite large in the initial down- 
stream region. From extensive4xperimentation, a smaller oonstant (K = 0.01 , equation (44)) was found 
to be effective for all three dynamic pressure ratios studied. Vhe turbulent Schmidt number distribution 
coincided with that of the single- and multiple-jet studies discussed. 


Shown in figuresf27 through 29 are the results of the virtual source computational simulation of 
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the multiple injection configuration of ref. 6 for dynamic pressure ratios (q r ) of 0.5, 1.0, and 1.5, re- 
spectively. Superimposed are appropriate single-jet data from ref. 5)\as well as the multiple-jet data. 

The computational results are compared, in the region of 0 < x/D < 60, on the.multiple bases of: 
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(a) maximum concentration level of hydrogen at x/D 
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(b) trajectory of the maximum hydrogen concentration level (i.e., 'transverse displacement 

/ above the plate of the local maximum hydrogen level) \ 

lateral spreading of the jet as determined by the z/D coordinate of\the hydrogen concentra- 
/ tion equal to 10% of the local maximum level. 
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Figure 28 COMOC Predictions for Virtual Source. q r = 1 

Scramjet combustor systems that may not be predominantly diffusion controlled are cer- 
tainly conceivable. For practical design studies it then becomes necessary to perform computations 
for finite rate or equilibrium controlled combustion systems to establish a data base, as well as to 
validate the usefulness of mixing efficiency computations from cold flow studies. The COMOC pro- 
gram system is readily extensible to analysis of combustion systems by addition of state-of-the-art 
kinetics packages, in the same direct manner that the turbulent transport subroutine was adapted 
from an operational two-dimensional Finite difference computer program. This program (NUMINT) 
was extensively utilized as well, for the accuracy and convergence determination of COMOC solutions, 
using the single-jet centerplane data of Rogers (ref. 5) as an initial condition (refer to Appendix A). 


38 








Longitudinal Displacement - x/D 


Figure 29. COMOC Predictions for Virtual Source, q f = 1.5 

To illustrate the near term potential of a practical combustion calculation capability, the temperature 
of these same initial data was elevated to autoignition, and the centerplane run repeated with finite 
rate combustion occurring. Shown in figure 30 are the computed hydrogen mass fraction distributions 
at x/D = 30 for the chemically reacting system and a comparison frozen flow computation as well as 
the initial distribution at x/D = 7. 
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Figure 30. Predicted Hydrogen Mass Fraction Distributions at x/D = 30 Using Experimental 
Centerplane Initial Data of Rogers (ref.5) 

These computations serve to illustrate the operational capability for numerical analysis in 
chemically reacting systems. Since subroutines are easily adapted for use by COMOC (as a conse- 
quence of its modular structure), a combustion calculation capability for practical three-dimensional 
flow fields could be readily established by borrowing from the state of the art. The prospect for 
practical design guidance of complex combustors thus appears achievable in the near term, and 
COMOC would benefit directly from accrued pertinent experience, in particular, the quasi-linear- 
tion approach of Morretti (ref. 30), used in conjunction with the implicit integration algorithm of 
Tyson and Kliegel (ref. 31). These important features are a part of the kinetics analysis capability 
(ref. 32) producing the example results. 
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CONCLUDING REMARKS 


A theoretical and analytical study has been made of mixing in the three-dimensional, turbu- 
lent boundary layer-type flow field downstream of transverse injection. Numerical solutions of the 
flow field equations were achieved using a general purpose computer program (COMOC) founded 
upon a finite element solution algorithm. Computations using a turbulent transport model based on 
mass defect and mixing length concepts, in conjunction with a geometrically defined turbulent Schmidt 
number, yielded reasonable comparisons with experimental data for both single and multiple injector 
geometric configurations. Refinement of the model is certainly required before any generality may 
be assigned to it. 

Detailed characterization of the complex three-dimensional flow field has resulted from 
this analytical study. The concept of a virtual source was evaluated successfully, and mixing effi- 
ciency computations were illustrated as potentially useful for support of engineering design studies. 

Of major consequence may be the observation that agreement with data occurred only for the 
locally frozen hydrogen concentration boundary condition directly beneath the jet centroid. This 
isolation from the wall may well exert great influence on combustor design and operation. 

It is clear that the controlling physics of such phenomena cannot be established by numerical 
observation alone. However, the COMOC computer program system, with its demonstrated flexi- 
bility and versatility, may have considerable impact on the bringing together of the present and next 
generation of digital computers with the engineer and designer. The theoretician will certainly re- 
ceive guidance in his work using a versatile computational capability as his laboratory, whereas the 
experimentalist can employ the computer to ascertain flow regions requiring especially detailed 
measurements by rapid comparison of data and computations. 
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APPENDIX A 

ACCURACY AND CONVERGENCE OF THE FINITE ELEMENT 
SOLUTION ALGORITHM IN COMOC 


Favorable evaluation of numerical accuracy, and convergence with discretization, of 
a COMOC-generated finite element solution for the two-dimensional steady Navier-Stokes 
equations is reported (ref. 15). Accuracy and convergence for the transient intergration of 
the finite element solution algorithm for heat conduction is proven (ref. 13). This Appendix 
discusses results which establish a favorable assessment of COMOC-generated solution accuracy 
and convergence for boundary layer flow of a compressible fluid. 

Comparison solutions were generated using an operational finite difference computer 
program (NUMINT) which integrates the equations governing multi-component mixing and 
reaction of compressible, two-dimensional and axisymmetric boundary layer flow, cast in the 
von Mises coordinate system. The use of derived variables in NUMINT, as distinct from the 
physical variables of COMOC, renders available two completely independent methods for 
solution of the governing equation system. NUMINT has been extensively evaluated, and is 
in everyday use for prediction of practical mixing and reacting compressible boundary layer 
flow fields. 

The symmetry plane, parallel to the freestream velocity vector and bisecting the 
injection orifice, figure 2 was selected as the sample comparison solution domain, since NUMINT 
is constrained to two-dimensional problems. The three-dimensional capability of COMOC was 
numerically reduced to two by enforcing vanishing gradient boundary conditions on all 
dependent variables on two planes, z = ± constant, and prescribing initial conditions that were 
independent of z, i.e., uniform across the discretization. By this procedure, COMOC was 
independently assessed to preserve two-dimensionality within about ±2% by application only 
of this non-rigid boundary condition statement. 

The single jet centerplane data of Rogers (ref. 5), for q r - 1 .0 and at Station x/D = 7, 
was selected as the initial condition for longitudinal velocity, u, and hydrogen mass fraction, 

Y. The cross flow velocity, w, was set uniformly to zero, and COMOC started computations 
for transverse mass flow pv after sufficient data were generated to allow approximation of 
(pu( , equation (31). Vanishing normal gradient boundary conditions were applied on the 
closure of the elliptic domain except for enforcing no-slip at the plate surface. Computations 
were carried to Station x/D = 30, and to stimulate strong mixing and to control parameter 
variation, a constant turbulent eddy mass mixing coefficient was specified, and the Schmidt 
Number was set constant throughout, equal to 0.7. 

Shown in figure A.l are comparison COMOC and NUMINT computed solutions for 
hydrogen mass fraction distribution through the boundary layer thickness at Station x/D = 30. 
Extensive diffusion and convection has occurred, as illustrated by comparison to the initial 
distribution (x/D = 7). A one-to-one-correspondence exists between COMOC data points 
and the discretization through the boundary layer thickness. The curve is faired through the 
square COMOC symbols corresponding to the 22 finite element discretization. Agreement 
with the 100 zone finite difference NUMINT solution is excellent, with only the values of Y 
near the wall being computed as measureably different. NUMINT employs a four-point finite 
difference derivative formula to evaluate Y at the wall to enforce the vanishing gradient 
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boundary condition. COMOC, of course, handles Y at the wall in a fashion identical to interior 
node points, resulting in simultaneous evaluation of the wall value of Y with the interior 
solution. The small computed differences might well stem from the distinctly different 
methods for handling such boundary conditions. 

The important proof of convergence with discretization is also illustrated in figure A.l . 
The open circles correspond to the COMOC solution using a very coarse 13-finite element 
discretization. Agreement with the 22-element solution is amazingly good, especially when 
considering that the initial peak Y distribution at Station x/D = 7 was captured at only one 
node. More importantly, the 13-element solution is observed to depart from the 22-element 
solution in directions diametrically opposite to that of the 100-zone finite difference solution. 

A 50-zone finite difference solution, also generated by NUMINT, agreed almost exactly with 
the 100-zone solution, except at the plate surface, where the inverted triangle data points 
confirm the identical solution diametric departure trend observed. From the mathematics 
standpoint, since the finite difference solution algorithm is amply proven to converge with 
discretization, the same property can be ascribed to the finite element algorithm. From the 
practical analysis viewpoint, the excellent solution accuracy for the 13-finite element discreti- 
zation indicates the ability to employ rather coarse discretizations, and yet obtain adequately 
accurate numerical solutions for practical problems. Additional detail on this subject is 
presented in reference 13. 

The comparison of the computed longitudinal velocity distributions is shown in figure 
A. 2. A strong retardation of the flow is indicated, by comparison to the initial velocity 
distribution, and negligible differences exist between the COMOC and NUMINT solutions for 
22-finite elements and 100-zone finite difference, respectively. The COMOC solution for the 
13-element discretization was essentially identical to the 22-element case. 

Two versions of a global continuity equation solver are operational in COMOC. The 
finite element approximation produces polynomials in the y coordinate only, and analytic 
integration of the solution algorithm, equation (30), is directly obtained. The distinctive 
feature of the two algorithms is the finite element approximation to the longitudinal mass 
flux dependent variable derivative vector |pu[ , equation (31 ). Both running-smoothing 
polynomial representations, over sequential panels of data, and cubic spline fits over the entire 
data field have yielded excellent results. Shown in figure A. 3 is the comparison between COMOC- 
generated solutions using quadratic running-smoothing (over three data points) and the cubic 
spline, for transverse mass flux distribution, pV, obtained for the axisymmetric companion 
problem to center plane mixing parallel hydrogen injection into a supersonic air stream. The 
identical radial distribution of JpUj was used for each case, and agreement between the 
solutions is excellent. Slight underprediction by NUMINT, in the back flow region, just off 
the flow centerline is indicated. This may result from large truncation errors associated with 
numerical differentiation of streamline data over a fine discretization near the flow centerline. 

An important and independent assessment of computed solution accuracy is numerical 
evaluation of properties of the flow field that would be rigorously conserved if an analytic 
solution to the governing differential equation system could be achieved. As an example, the 
species continuity equation for hydrogen mass fraction, equation (4), can be written in 
explicit conservation form, using the three-dimensional gradient operator v, as 

V (pUY-^ 7 2 Y) = 0 (A.l) 
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APPENDIX A 


\ACCURACY AND CONVERGENCE OF THE FINITE ELEMENT 
SOLUTION ALGORITHM IN COMOC 



Favorable evaluation of numerical accuracy, and convergence vyith discretization, of 
a COMOC-generated finite element solution for the two-dimensional^steady Navier-Stokes 
equations is reported (ref. 15). Accuracy and convergence for the transient intergration of 
the finite elemenbysolution algorithm for heat conduction is proven 7 (ref. 1 3). This Appendix 
discusses results which establish a favorable assessment of COMOC-generated solution accuracy 
and convergence for\>oundary layer flow of a compressible fluid. 

Comparison solutions were generated using an operational finite difference computer 
program (NUMINT) which integrates the equations governing multi-component mixing and 
reaction of compressible, two-dimensional and axisymmetric boundary layer flow, cast in the 
von Mises coordinate systenk The use of derived variable's in NUMINT, as distinct from the 
physical variables of COMOcXrenders available two completely independent methods for 
solution of the governing equation system. NUMINT/has been extensively evaluated, and is 
in everyday use for prediction ofxpractical mixing and reacting compressible boundary layer 
flow fields. \ / 


The symmetry plane, parallel to the freestream velocity vector and bisecting the 
injection orifice, figure 2 was selected asSthe sample comparison solution domain, since NUMINT 
is constrained to two-dimensional probl^ms.yThe three-dimensional capability of COMOC was 
numerically reduced to two by enforcing vanishing gradient boundary conditions on all 
dependent variables on two planes, z = ± constant, and prescribing initial conditions that were 
independent of z, i.e., uniform across the'discretization. By this procedure, COMOC was 
independently assessed to preserve tw^rclimen«onality within about ±2% by application only 
of this non-rigid boundary condition statement.^ 

The single jet centerplane data of Rogers (ref. 5), for q r = 1.0 and at Station x/D = 7, 
was selected as the initial condition for longitudinal velocity, u, and hydrogen mass fraction, 

Y. The cross flow velocity, w, was set uniformly to zeto, and COMOC started computations 
for transverse mass flow pv aft^ sufficient data were generated to allow approximation of 
|pu} , equation (31). Vanishing normal gradient boundary conditions were applied on the 
closure of the elliptic domaim except for enforcing no-slip Vt the plate surface. Computations 
were carried to Station x/dA 30, and to stimulate strong mixing and to control parameter 
variation, a constant turbulent eddy mass mixing coefficient \as specified, and the Schmidt 
Number was set constantXhroughout, equal to 0.7. 

Shown in figure A.l are comparison COMOC and NUMINT computed solutions for 
hydrogen mass fraction distribution through the boundary layer thickness at Station x/D = 30. 
Extensive diffusion /nd convection has occurred, as illustrated by comparison to the initial 
distribution (x/D ^ 7). A one-to-one-correspondence exists between QOMOC data points 
and the discretization through the boundary layer thickness. The curvets faired through the 
square COMOC/symbols corresponding to the 22 finite element discretization. Agreement 
with the 100 Zone finite difference NUMINT solution is excellent, with only the values of Y 
near the waj/being computed as measureably different. NUMINT employ^a four-point finite 
difference derivative formula to evaluate Y at the wall to enforce the vanishing gradient 
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Hydrogen Mass Fraction - V (%) 



Vertical Displacement - y/D 


Figure A.l . Comparison Computed Hydrogen Mass Fraction Distributions on Symmetry Centerplane 

of Mixing Region 


Identify a three-dimensional, rectangular parallelepiped control volume, figure A.4, with sides 
normal to coordinate directions. Employ Gauss’ theorem to recast equation (A.l) as a surface 
integral over this control volume: 
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Figure A. 2. Comparison Computed Longitudinal Velocity Distributions on Symmetry 

Centerplane of Mixing Region 


Integrand evaluation over tlie subject control volume yields, alternatively. 
puYda ~ P l 'Yda + pvYda 


k 


ii Le 

Y. da + 
Pr z 


fD 


pLe r u Le 

- — Y, v d a+ - — Y. da 
Pr y I Pr y 


- 0 


(A. 3) 


In equation (A. 3), the integral subscripts identify control volume surfaces, and the symmetry 
properties on face E and the no slip wall have been accounted for. 
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Hydrogen Mass Fraction - V (%) 
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f Vertical Displacement^y/D 

Figure A.l . Comparison Computet! Hydrogen Mass Fraction Distributions on Symmetry Centerplane 

/ of Mixing Region \ 


Identify a three-dimensional, rectangular parallelepiped control volume, figure A.4, with sides 
normal to coordinate directions. Employ Gauss’ theorem to recast equation (A.l) as a surface 
integral over this control volume: \ 


(A. 2) 


)[pUY - v 2 Y] • nda = 0 






Transverse Mass Flux (g/Lt) 


kg/m 2 sec 
25.0 r 


20.0 L 


15.0 


10.0 


5.0 


0 U 


-5.0 L 


- 10.0 


Ibm/ft 2 sec 



0.002 0.004 0.006 0.008 0.010 0.012 0.014 0.016 0.018 0.020 (m) 


Radial Displacement (r) 

Figure A. 3. Computed Transverse Mass Flux Distribution Comparison for COMOC and 

NUMINT Solutions 

For all computations reported herein, the hydrogen mass fraction vanishes on face D. 
Selective terms may vanish from within the square bracket, equation (A. 3), dependent upon 
applied boundary conditions, e.g., Y, z vanishes on surface C for the multijet, virtual source 
and centerplane computations, while Y,y essentially vanishes on surface D for all computa- 
tions, since the solution domain extends well into the freestream. 
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Figure A-4. Control Volume for Numerical Evaluation of Conservation Properties of 
Three-Dimensional Flow Fields Computations by COMOC. 

If vanishing normal gradient boundary conditions are enforced at the plate surface 
as well, the conservation expression reduces to the familiar form 

/ puYda = constant (A. 4) 

J A 

which states that the hydrogen mass fraction flux is analytically constant in every plane witli 
normal parallel to the predominant flow direction. The approximate adherence of computed 
solutions to equation (A. 4) is considered important for verification of the proper operation 
of a computer program. 


The finite element solution algorithm provides a mathematically consistent and 
readily obtained numerical approximation to equation (A. 4). From equation (21), written 
on both pu and Y. obtain the following approximation to equation (A. 4) within a finite 


element. 

FLOW m 



|x||x| T do [ T 1 ni (Y) m 


(A. 5) 


The integral in equation (A. 5) is readily evaluated. A simple DO loop over the finite elements 
of the discretization sums M terms of this form to yield a numerical approximation to 
equation (A. 4) which contains the details of the computed distribution of both pu and Y. 
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Figure A. 3. Computed Transverse Mass Flux Distribution Comparison for COMOC and 
/ NUMINT Solutions ^ 

For all computations reported herein, the hydrogen mass fraction vanishes on face D. 

Selective terms may vanish from within the square bracket, equation^ A. 3), dependent upon 
applied boundary conditions, e.g., Y, z vanishes on surface C for the multijet, virtual source 
and cente^ane computations, while Y,y essentially vanishes on surfaceT) for all computa- 
tions, sin^e the solution domain extends well into the freestream. \ 
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COMOC computes the sum of equation (A. 5) at every output station. Excluding a gross 
programming error or inconsistent theory, discretization exerts the predominant influence 
on the sum. In actuality, solution surfaces are complexly curved; equation (A. 5) performs a 
piecewise planar interpolation of the surface, yielding a multi-faceted approximation to the 
true surface. 

The previously discussed centerplane computations employ the correct boundary 
condition distribution for equation (A. 4) to be rigorously valid. Shown in figure A.5 are 
COMOC-computed deviations from absolute hydrogen flow conservation for the 13- and 22- 
finite element solutions of the centerplane problem. The comparison data base in both tests 
was total computed hydrogen flow at Station x/D = 7. The 22-element COMOC solution 
numerically conserves total hydrogen flow to within ± 1.5 percent. There is a computed 
continual loss of hydrogen mass flow for the 13-element solution, amounting to about 10 
percent of the original computed flow at Station x/D = 30. Referring to figure A.l , this 
computed hydrogen flow loss corresponds to a maximum solution deviation for hydrogen 
mass fraction of about 3 percent, occurring at the peak of the initial input distribution, to 
about -4 percent, occurring at the plate surface, and near freestream. Considerable additional 
experimentation is, of source, required to establish a firm correspondence between computed 
hydrogen flow conservation and actual solution accuracy. At this juncture, a 10 percent 
computed loss appears as an acceptable indication of adequate solution accuracy. 



for Different Discretizations 


The centerplane tests were repeated with hydrogen mass fraction held constant at the 
initial value at the wall, for 7 < x/D < 30. The maximum computed deviation from hydrogen 
conservation was +0, -2 percent. Thus, equation (A. 4) is probably reasonably accurate for 
those cases in the present analysis where the diffusion terms (in square brackets), equation 
(A. 3), do not rigorously vanish by gradient boundary condition statements. 
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APPENDIX B 


VIRTUAL SOURCE INITIAL CONDITIONS 


The complexity of flow in the near injection region requires that detailed initial data distri- 
butions be known to start computations at a downstream station where the elliptic boundary layer 
approximation is valid. The alternative is establishment of an analytical and/or empirical model of a 
numerical “virtual source”, as discussed. This approach requires establishing simplified initial condi- 
tions at a station upstream of the mixing region as a function (only) of undisturbed freestream and in- 
jection parameters. 

Injection of a jet from an orifice in a plate into a transverse supersonic air stream has been the 
subject of a number of investigations (ref. 25 through 28). An important correlating parameter is q r , de- 
fined as the ratio of the dynamic pressure in the jet to the corresponding freestream value. Except for 
the investigations by Rogers (ref. 5 and 6), available experimental data are typically for large values of 
q r , whereby the jet has sufficient momentum to penetrate the boundary layer and produce the com- : 
plicated separation region and bow shock ahead of the jet. In references 5 and 6, q r ranges between 
0.5 and 1.5. The jet lacks the necessary momentum to penetrate the boundary layer. Therefore, mix- 
ing occurs within a turbulent boundary layer velocity profile. Consequently, the referenced empirical 
models are not directly applicable to analysis of the present data, and other means of characterizing 
the near injection region were investigated. 

The proposed barrel shock model of the turning jet is shown in figure B.l . For large values of 
q r , a similar configuration with an interaction bow shock- has been considered in reference 25. An 
analysis based on one-dimensional flow was developed for the present case of small q r to characterize 
the turning jet. The parameters for the present model are listed in table B.l. 



, < 2 > 

Figure B. 1 . Transverse Injection into a Turbulent Boundary Layer 
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TABLE B.l 

BARREL SHOCK PARAMETERS 


7 

specific heat ratio 

W A = 

molecular weight of freestream (air) 

m 2 = 

local Mach number of region 2 

w H = 

molecular weight of jet (hydrogen) 

Pti = 

total pressure in region 1 

C 1 = 

orifice discharge coefficient 

Pt3 = 

total pressure in region 3 

w 4 = 

orifice diameter 

Ttl = 
T t3 = 

total temperature in region 1 
total temperature in region 3 

m 5 = 

local Mach number in region 5 


The effective diameter of the stream tube that passes through the Mach disk is described by a dis- 
charge coefficient. The isentropic one-dimensional flow relations, and the condition of a sonic jet 
throat, M 4 = 1 , are used to determine conditions at stations 3 through 5. The normal shock relations 
yield the flow quantities in region 6 , and the jet is expanded in region 7. 

The analysis procedure is an iteration to find the Mach number upstream of the Mach disk, 
M 5 , that allows jet static pressure equilibration with local freestream, i.e., p 7 = P 2 - From M 5 and the 
normal shock relations, the velocity immediately downstream of the Mach disk is determined (Ug). 
Illustrated in figure B.2is the determination of Ug for Roger’s test conditions of q r = 0.5 and q r = 1.0; 
assuming Cj = 1 .0, i.e., all the injectant passes through the Mach disk. For the virtual source compu- 
tations discussed, the initial velocity within the hydrogen jet was identified with U^, as determined by 
this procedure. 

A FORTRAN listing of the barrel shock computer program is included as figure B.3. The re- 
quired input for a problem is shown in table B.2, and a representative output sample is attached as 
figure B.4. 


TABLE B.2 

INPUT SPECIFICATIONS FOR BARREL SHOCK COMPUTER PROGRAM 


Card No. 

FORTRAN Input Variables* 

Format 

1 

GAMMA, M2, PTI 

3F10.5 

2 

PT3, TT1 , TT3 

3F10.5 

3 

WA, WH,C1 

3F10.5 

4 

W4 

FI 0.5 

5 

XM5 

F10.5 

6 

NCASE 

13 


*See Table B. 1 
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APPENDIX B 


VIRTUAL SOURCE INITIAL CONDITIONS 


\ The complexity of flow in the near injection region requires that detailedinitial data distri- 
butions be known to start computations at a downstream station where the elliptic boundary layer 
approximation is valid. The alternative is establishment of an analytical and/£r empirical model of a 
numerical “virtual source”, as discussed. This approach requires establishing simplified initial condi- 
tions at a station upstream of the mixing region as a function (only) of u/ldisturbed freestream and in- 
jection parameters. 
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lnjectio^i of a jet from an orifice in a plate into a transverseyupersonic air stream has been the 
subject of a number of investigations (ref. 25 through 28). An important correlating parameter is q r , de- 
fined as the ratio^the dynamic pressure in the jet to the corresponding freestream value. Except for 
the investigations by, Rogers (ref. 5 and 6), available experimental data are typically for large values of 
q r , whereby the jet has sufficient momentum to penetrate the4>oundary layer and produce the com- 
plicated separation regi’on and bow shock ahead of the jet. In references 5 and 6, q r ranges between 
0.5 and 1.5. The jet lacks the necessary momentum to penetrate the boundary layer. Therefore, mix- 
ing occurs within a turbulent boundary layer velocity profile. Consequently, the referenced empirical 
models are not directly applicable to analysis of the present data, and other means of characterizing 
the near injection region wer^investigated. 
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The proposed barrel shock^model of the turning jet is shown in figure B.l . For large values of 
q r , a similar configuration with an interaction bovyfhockhas been considered in reference 25. An 
analysis based on one-dimensional flo,w was developed for the present case of small q r to characterize 
the turning jet. The parameters for the^present model are listed in table B.l. 
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Figure B. 1 . Transverse Injection into a Turbulent Boundary Layer 
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Velocity U 








Figure B.3. Listing of Barrel Shock Computer Program 

HEAD (5,101) G , XM 2 , PTTTpT 3 ', TT 1 , TT 3' , V7A , T-7} T , C 1 , W 4 

22 WRITE (6,204) • - ■ ■ 

20 4 FORMAT (1111 ) 

READ (5,101 ) XM5 

101 FORMAT ( 3F 1 0 . 5 ) 

_ P-r’ AD ( 5_,_16 ) MCASX 

16 FORMAT (13) 

G 1 = G + 1 . 

~*G2 = G - T. 

G3 = 0.5 * G1 
G4 = 0.5* 02 

G5 =__2_. * _G/G1 


XH5S 
EX 7 = 

= XM5 * XII 5 
G3/G2 


EX o = 

-EX 7 


AS1 = 

G3 ** EX 7 


AS 2 = 

1 . + G4 * XM5S 


AS 3 = 

AS 2 ** EX0 


AS 4 = 

AS 1 * XII 5 * AS 3 


S50S4P = 1./AS4 


XM2S 

= XI! 2 * X.M2 


RA = 

1545. /VIA 


RII = 

1545./WH 


XX = 

RA * TT 1 


DTI = 

(144. * PT1)/XX 


DT3 = 

(144. * PT3 ) / (RM 

* TT 3 ) 

7,2 = 

1 . + G4 * XM2S 


El = 

-1 . 


E2 = 

-G/G2 ' 


E3 = 

-1 ./Cl 2 


E4 = 

-0.5 


XK1 = 

7 2 * * E 1 


II 

72**E2 


XK3 = 

72 ** E3 


XK4 = 

72 ** E4 


TT2 = 

TT 1 


PT2 = 

TT1 


DT2 = 

DTI 


AT IS = 

= G * 32.2 * RA * 

TT 1 

ATI = 

SQPT (ATI S) 


AT 2 = 

AT 1 


AT3S = 

= G * 32.2 * RII * 

miyi O 
j. _ — » 


AT 3 = SQRT ( AT 3 0 ) 
T2 = TT2 * XK1 
P2 = PT2 * XK2 
D2 = DT2 * XK3 


A 2 

= AT 2 * 

XK4 


V2 

= AT 2 * 

XM2/(SQRT (72) ) 


TT4 

= TT3 



PT4 

= PT3 




DT4 = DT3 
AT 4 = AT 3 
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Figure B.3. —Continued j 

XM4 =1. 

7,4 = 1 . + G4 * XM4 * XM4 
Q.1 '="7 4 **E1“ 

Q2 = 74 ** E2 , 

-~Q3~= -ZTTWK3 s 

Q 4 = 74 ** E4 
T4 = TT4 * Q1 

P4 = _PT4_* Q2 

D 4 ~= DT4 * Q3 ““ 

,A4 = AT 4 * Q4 ' ^ V 

V 4 = AT 4 ♦ XM4/TSQRT (74) ) 

QR = D4_ * V4 * V4/(D 2 * V2 * V2) _____ 

PI = 3 7 1 416 ■ . f V ‘ ■ v T-' -■ 

JI4 1 = PI _*_W4 * W4 /(4 . * 14 4.) _ 

F4~= D4 * V4*~'* _ H4 

F4il =453 ,_6 *__F4 

EXT = G/G2 

EX2 = J./G2 ' . v'y ■' ; . . . _ • 

W4P = Cl * W4 

CS4 = 4_./PI _ 

S4P = W4P * W4P / CS4 

S5 = S4P * S50S4F _ __ ; ? ___ ________ 

W5S =' CS4 * S5 ' - 

W5 .= SQRT(W5S)_ 

XH5S~=XM5 * XMff . ' " “ 

75 = 1 G4 * XI I 5 S . . . 

' AQ 1 " = 7 5 * * E 1 

AQ 2 = Z 5 ?_*_ E2 . 

AQ3 = 75 ** E3 ' ' ~ ” ~ 

AQ4 = 75 ** E4 ;* , ^ . : 

~TT 5 = TT4 * ” 

PT 5 = j?T4__ 

DT5 = DT4 " ‘ . * “ 

_AT5 = AT 4 _ 

T5 = TT5~AQ1 

P5 = PT5 * 7 \Q 2 ‘ 

D5 = DT5 * AQ.3 
A5 = AT 5 * AQ4 
*V5 = AT5 * XM5/ (SQRT (75)T 

P60P5 =_ ( 2. * G * XH5S - G2)/G1 

D60D5 = G1 * XM5S/ (G2 * XM5S + 2.) "" ~ ‘ 

T60T5 = P60P5 * (G2 * XM5S + 2.)/(Gl * XM5S) 

V60V5 = 1 . -'4. * (XM5S - 1.) * (G * XM5S +" 1 .T7TWI **XH5S')^ r * 2 .7 

A 60 A 5 = SQRT (T60T 51 ' __ . 

' AS V <33 * XT15S * - 

BS = 1./P 60P5 

TT6 = TT5 ’ ” 

JPT6_ = P5 * (AS ** EX 1 ) * (BS ** EX2 ) 

DT6 =~ PT6 * 14 4.’ * WH/(1545. * TT6 ) ’ ’ 

AT 6 = AT 5 ; . . _ 

~ W6‘ = W5 
T6 = T5 * T60T5 

P6 = P5 * P60P5 ‘ 

D6 = L)5 _* D6QD5 
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Figure B.3. -Continued 

A 6 = A5 * AGOA5 
DEEP = P6 - P2 

I F ( DELP . LT .0.0) "GO TlT'l'TS 1 ' : 

XM6S = (G2 * XM 5 J i + 2.)/ (2. * G * XM5S - G2) 

"' XIIGT = SQRT (XM6S) ' 

7C = 1. + G4 * XI! G 5 

vrr = AT 6 * XKB/'(SQRTTfT(5T) 

TT7 = TT6 

PT7 = PT6 

DT7 = UTG 

AT 7 = AT 6 - • - - 

P 7- = P2 

PF = P7/PT7 

PFG = PF ** (-G2/G) 

XM7S = (PFG - T.) /G4 ’ " ’ 7 

JIM 7 = SQRT (XI17S ) 

XH7S = XM7'*'~Xn7' ‘ 

7-.1 = 1 . + G4 * XI! 7 S 
UQ 1 = 7.1 * *E1 

BQ2 = 71 _** E2 

BO 3 = 71 ** E3 
BQ 4 = 7,1 ** E4 
T7 = TT7 * BQ1 

P7C = PT7 * BQ2 _ 

D7 = DT7 * BQ 3 

A 7 = AT 7 * BO 4 

V7 = AT 7 * XI 17/ (FQRT ( 7,7 ) ) 

CJI; = 1 . + G 4 * XH6S 
CJD = l7' + G4 * Xf-I7S ‘ 

CJE = CJN/CJI) 

CJF = III16/XH7 

CJG = -G3/G2 

CJII = CJF * (CJE ** CJG) 

CJI = SQRT (CJII) 

V; 7 = WC * CJI 
FCOI1 = PI/576 . 

F4 = U4 * V 4 * 174 '* >74 * FCON 

F 4P = U 4 * V4 * 17 4 P * II 4 P * FCOE 

F5 = 05 * V5 * 05 * 175 * FCON 

FG = DG * VG, * 176 * 176 * FCON 

F7 = 07 '* V7 * 177 * IJ7 * FCON 

1 1 8 CONTINUE 

1 'RITE (6,12 0 ) NC/vSE 

120 FORMAT (T10 , ' SECONDARY JET ' , / , T 1 0 , ' prOGRA?!- TP PY F.O.IIAINS A 

* . OKULEV7ICn ' ,//// , T 1 0 , ' CASE NO . ' 13) 

WRITE (6,121) G,XH2 

121 FORMAT (////,T1 0 , ' INPUT ',///, TIP , 'G' ,T 1 ° , r ' = ' ,E12.U,T42, 

* ' SPECIFIC HEAT RATIO' , /TIP, 'XM2' , T 1_ n ,7_=7_,_F 1_2^4 ,T4 2 , 

* ' I lACII NO.’ IN “REGION 2 ,V REE STREAM r ) •’ 

WRITE (6,122) _PT>, PT3 ,_TT 1 

122 FORMAT (T10, r PTl ' ,T19, ' = 1 , El 2 . 4 ,T42 , ' TOTAL PRESSURE , P-PGiOT- 1 
*T82 , ' (rSIA) » , / , T1 0 , ' PT3 • ,T1°, ' = ’ ,E12.4,T42, ' TOTAL PPESST^p^ 




p 5 = PT5 *_ AO 2_ / \ 

55 = DT5 ♦ AQ3 "/ 

A5 = AT 5 * AQ4 / 

V5 = AT5 XM5//( SQRT ( 7 , 5 )T~ 

p 6 OP 5 = ( 2. */G * XM5S - G2)/G1 

D60D5 = G1 * XM5S/(G2 * XM5S + 2.) 

T60T5 = P60P5'' * (G2 * XM5S + 2.)/(G1 

V60V5 = 1 . -f 4.'* (XM5S - 1.) * (G * XM5S\+ “1 .') /T(GY* *XM5"S ) '*'* 2 .7 
A60A5 = SORT (T6QT 5) 

AS ="'Gy“+XH5S 

IiS_= 1 . /P 6'O P 5 

VfiT< g = >jirr Cy> 

PT6 = P5/* (AS * * EX1 ) * (BS ** EX2) 

DT6 * 144.* WH/ (1545. * TT6 ) 

AT 6 = /AT 5 , . 

W6"= V75 ' 

T 6 = T5 * T60T5 
P6 = P5 * P60P5 

1)6 = D5 _* D 60D5_ 

'' '55 





Figure B. 3. — Continued i. 

' ’ 

*T7t'827 t TpSIA) • ,/,Tl 0 ,'TT1 ' ,'.T 1 9 = ' ,E12 1 ; 4/,T42 , ’.TOTAL 7 TEMPEPATURjf,TuT 

* GljOtJ 1 1 , T 8 2 , 1 (DEG . RANKINE) ') > _ '*.«*< L5>:; \ j . __ __ 

WRITE (6, 123) TT 3 , WA , WII 

123 FORMAT (T1 0, 'TT3' # Tl 9, ' = ' ,E12.4,T42, 'TOTAL TEMPERATURE REGION 3', 

* " *T 8 2 , (DEG . "RANK IliE ).',/, T 1 0'^A ' ,.T 1 9 , ' = ' ,E1 2 . 4,T42, 'MOLECULAR" WEIGH"' 
* T/ A1R' f T82jt ' ( 1/M OLE) 1 , /VTl 0 , * WH ' ;T19 ;-&»«■£ El 2 ffr,T42r*MOEECULAR; ; W EIG 

*IIT, HYDROGEN' ,T82 # 1 (1/MOLE) ') » 

WRITE (6,. 124) Cl ,W 4,XM5 

12 4" FORMAT (T 1 0 , * C 1 ' , T1 9. , ' = ' ,E1 2 . 4 ,T42 , ' ORIFICE DISCHARGE COEFFICIENT ’ 
*J__ " L. ; . 

♦ , T10 , ' W4 ' ,Tl 9 , ' = ' ,El 2 . 4 ,T42 , 1 DIAMETER OF ORIFICE ' ,T8?,' , llNCHES)“ r ; 
*/,Tl 0 , 'XM5 ' ,T1 9 , ' = ' , F.12.4,T42, 'MACH NO. IN REGION 5') 

IF(DELP.LT.OVO) GO TO 202; 

WRITE (6 . 125) PT1 , DTl ,TT1 : v' : 

1 2 5 '"FORMAT T/AtIO, 'OUTPUT '7//,Tl0, *PT1 ' ,T19 , ' = ' ,E12. 4,T42~ 

* 'TOTAL PRESS 

"*'URE , REGION 1 ' ,T82,' (PSIA) ' ,/,TlO', • DTl^yTl?, * = ' , El 2 . 4 TOTAL DE" 

♦NSITY , REG I OIL 1 ' ,T 82 , ' (LBM/ CU . FT . ) ' , A T1 0 1 TT1 ' , T 1 9 , ' = 1 ,E12 .4,T42 , 

* 'TOTAL TEMPERATURE , REGION 1 ' ,T82 , ' (DEG. RANKINE) ' ) 

WRITE (6, 126) ATI ,PT2 ,DT2 

126 FORMAT (tYO, ' AT 1 "' , T 1 9 , ' = ' , E 1 2 . 4 , T 4 2 , ' TOTAL ’SPEED OF SOUND , REGION T' 

* , T 82 , ' (FT. /SEC.) ' ,AT 1 0; , PT2' ,T19, ' , E 1 2 . 4 , T 4.2 ' TOT AL PRESSURE , RE 

♦GION 2 '7T82~, ' (PSIA) ' ,"/’» T10 , 'DT2 ' ,T19 , '= ' , El 2 . 4 ,T42 , 'TOTAL DENSITY," 

* REGION 2 ' ,T8 2 , ' (LBM/CU . FT ♦ ) ' ) 

WRITE (6 127) TT2 AT 2 P2 

127 FORMAT (TIP.' TT2 ' \ T 1 9 [ ' = ' yE 1 2 . 4 , T 4 2 , ' TOTAL TEMPERATURE REG TP TT 2 ', 

*T82, ' (DEG . 'RANK I HE ) ' , /, T 10,' AT2 ' , T 1 9 , ' = ' ,E T2 . 4 ,T42 , 'TOTAL' SPEED OF ” 

* SpUIID , REG ION 2 ' , 182, ' (FT. /SEC.) ' , /,T1 0 , ' P2 ' ;T l 9 , , E 12.4, T42 f 

♦ 'PRESSURE' REGION 2 ' , T 8 2 , ' (PSIA) ' ) ' . 

WRITE (6, 128) D2 ,T2 ,A2 

1Z8 FORflATITT07’ r D'7 t 7TTT7 T '= T ,'E'T2' . 47T4T, ' DFNS IT v VRFGIOM~7” r 7 r ! 1 877 ' (LBM7CU7'~ 
♦FT.) ' ,/,TlO, ' T2 ' ,T1 9 , ' = ' ,El 2 . 4 ,T42 , ' TEMPERATURE , REGION 2' ,T82, 

* ' (DEG. RANKINE) ' ,/,TlO, 'A2' ,T19, '=' ,E12 . 4,T42 , ' SPEED OF SOUND, REG T 

♦ON 2' ,T82 . 1 (F T. /S EC. ) ' ) 

WRITE (6, 129) V2,PT3,DT3 . 

1 2 9 FORMAT ( T 10 t ' V2 ' , T 19 , ' = ' ,-El2<4, T42,' VELOCITY . REGION 2' ,T8 2 , ' (FTL / S 
*EC .) ',/, T-10, ' PT3 ' ,T19, ' = ' , E12.4,T42,' TOTAL PRESSURE REGION 3','T82, : 

* ' (PSIA) ' ,/,TlO, ' DT3 ' ,T19 , '= ' ,E 12.4,T42, 'T OTAL DEN SITY , RE GION 3' , 

* T82, ' (LBM/CU . FT. ) ' ) ‘ ’ ' " ■ 

WRITE (6 . 1 30) TT3 . AT 3 , PT4 

130 FORMAT (T10 , 'TT3' ,T19 , '=' ,E12. 4 ,T42 , 'TOTAL TEMPERATURE , REGION 3', 
♦T82 ' (DEG. RANKINE ) ' ,/,TlO , ' AT 3 ' ,T 1 9 , ' = ' ,E12.4,T42, 'TOTAL SPE ED OF 

* SOUND , REGION 3 f ,T82, ' (FT. /SEC. ) '7/",TlO ,• PT4 ' ,T19 , ' = ' ,.E1 2 '. 4 ,T42 , " 

* ' TOTAL PRES SURE . REGION 4',T82. '(PS I A)') • ; _ 

WRITE (6 ,131) DT 4 , TT 4 , AT 4 

1 3 1 FORT1AT _ (TIP , ' DT 4 ' ,T19 A ' = ' ,E 12 . 4,T42 , ' TOTAL DEN S I T Y , R EG I O N 4 • , T82 , 

* ' (LBM/CU. FT. T' , / , T 1 0 , ' TT 4 ' , T 1 9 , '=' ,E12 . 4,T42 ; ' TOTAL TEMPERATURE ,R 
♦EGION 4' ,T82 , ' ( DEG. RANKINE) ' ,/,Tl0, 'AT 4.* ,T19 , ' = ' , E 12.4,T 42, 'TOTAL 

* SPEED OF SOUND, REGION 4' ,T82, ' (FT. /SEC.) ') 

WRITE (6. 132) P4.D4.T4 
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Figure B.3. - Continued 

132 FORMAT ( T1 0 , ' P 4 ') T 1 9 , ' ,E1 2 . 4 ,T42 , ' PRESSURE , REGION 4 ' f T8l7 l ”(PSIA j ' 

* / /.TIP . '04' .TIP, '=' ,E12 . 4,T42 , 'DENSITY, REGION 4* .782, [ (LBM./ CU .FT . 
* ) ' , T 1 0 , ' T4 ' , T 1 9 , ' = ' , E 1 2 . 4 , T4 2 , 1 TEMPERATURE , REGION 4', 

*j?_8 2 .' (DEG. RANKINE)') 


WRITE (6, 133) A4,V4,W4,W4P 

1 3 3 FORI IAT ( TJ 10 , 'A4* , T 1 ? . ' = ' ,E 1 2 . 4 ,T42 , ' SPEED O F ROUN D , RE G ION 4 • , T 8 2 , 

* ' (FT . /SEC.") ' , / , Tl'o", ' V4 ' , T1 9 , ' = • ,E 1 2 . 4 , T4 2 , 'VELOCITY , REGION 4 ' , T82 

* , ' (FT. /SEC . ) ' . / .TIP, * W 4 1 , T 1 9 , ' = ' , E 1 2 . 4 , T 4 2 , 'JET WI DTJj , RFC 1 4 ' , 

* T82 , ' ( INCHES )"'"",/ ,T1 0 , 'W4P ' ,T1°, '=' ,E12.4,T42, ' EFFECT iVE ' JET WIDTH 

*, REGION 4' ,T82, ' ( INCH ES) ' ) 


WRITE (6, 134) PT5,DT5,TT5 
1 3 4 FORMAT (TIP , 'PT5 ' ,719 , ' = ' ,E12 . 4,742 , ' TOTAL PRESSURE , REGION 5 ' , rn 89 
* ' (PSIA) '~ /,T1 0 , ' DT5 ' , T 1 9 , ' = ' , F. 1 2 .4 ,T42 , 'TOTAL DENSITY, RE CIO! 7 “ ' 


2 , ' ( LBM . /CU . FT .1 / , TIP# 'TT5 ' f Tl 9 ,_' = F 1 2 . 4 , T4 2 , ' TOTAL TPMPJ3RAT 

IJRE , REGION 5 ' ,7(52,' (DEG . RANKING) ' ) “ 

\ -RITE ( 6 . 1 3 5 ) AT 5 , P 5 , D 5 


135 FORMAT (T10, 'AT5' ,719, ' = ' ,E12. 4,742, 'TOTAL SPEED OF ROUND , REGION 
*,732, ' (FT. /EEC.) ' , /T1 0 , 1 P5 ' ,T 1 9 , ' = ' , F 1 2 . 4 , T4 2 , ' PRESSURE , REGION 7 


* TR2 , ' (PSIA)'' , / T 10, ' D5 ' , T1 9 , ' = ' ,E12. 4,742, ' DENSITY , REGION 


c * 


* TB2, ' (LBM./CU. FT.) ') 

7 ;EITE ( 6 , 1 3fj j T 5 , A5 , V5 

136 FORMAT (TIP, 'T5' ,T19, ' = ' , E 1 2 . 4 , T 4 2 , ' TEMPERATURE , REGION 5' ,TR2, 

*' (DEG.’ RANKINE) r ,/,T*1 0 , ' A5 ' ,T19 , ' = ' ,E 1 2 . 4 , T42 , ' SPEED OF ROUND , RFC I 
*01! 5 ' ,T82 , ' (FT. /SEC.) ' . / .T 1 P . ' V5 ' . 71 9 , ' = ' ,E1?.4,T42, ' VELOGITV , PEOI 
*011 5 ' ,T82 , ' (FT . /RFC . ) ' ) 

WRITE (6,137) _X!L5,W5 _ 

1 3 7 FOR] IAT (TIP,' XI 15 ' , T 1 0 , ' = ' , E 12. 4, T 42 , ' ’ IACE H*f *.BEP , REGION 5 ' , / , T 1 P , 

* 'US' ,719 , ' = ' ,E12 .>4,742, ' JET WIDTH , REGION 5 ' ,78? , ' (INCHES) ' ) 

WRITE (6,1 38 ) PTC , DT6 , 776 * 

1 3 B FORMAT (710,' PT_6 ' , TJ 0 , ' = ' , E 12 ._4 ,_T 4_2 ,/ TOTAL PRESSURE , REG ION 6 ’ , TP 3 , 

* ' (PSIA) V-/,Tio, ' DT 6 ' ,7197' = ' ,E12.“4 ,742 , 'TOTAL DEER IT v , REG 1 07 ’ 6 ' , 

* 782,' (LBM./CU. FT.) ' ,/, TlO, ' 776' ,T1R ’ = » ,rl?.4,T42 y 'TOTAL TFMPERAT 
*URK, REGION 6 ' ,732 , ' (DEG . RAMKINE) ' ) 

WRITE (6, 139)_A>T6,P6,D6 

139 FORILAT (TIP , 'ATG 1 ,719 ,"' ="* ,E12. 4,742, 'TOTAL SPEED OF ROUT’D , REGION 6' 
*,782, ' (FT. /SEC.) ' ,/TlO, 'P6' ,T19, '=' ,E12. 4,742, ' PRESSURE , REGION 6' , 

* 782,' (PSIA) ', /T1 0/' 1)6 ' jTlP , , = ' ,E12 . 4/T42, 'DENSITY /REGION 6', 

* 782, ' (LBM./CU. FT.) ' ) 

\ /RITE ( 6 , 1 4 P ) 76 , A 6 ', V6 

1 4 P FORMAT (TIP , '76 ' ,719 ,_' = ' ,E 1 2 . 4 ,T 42 , ' WE? IPERAT.UP.F , REGION 6 ' ,TR2, 

* ' (DEG. RANKINE) ' ,/,TlO ,"'A6*' ,719,'=' ,E12A,T42, 'SPEED OF 'SOUND*, PEG I 
*01! 6' ,TR2, ' (FT. /SEC.) ' ,/,TlO, 'V6' ,T1°, •=• ,F12. 4,742, ' VELOCITY , REG I 

* ON 6 ' , T 82 , ' ( FT . /SEC", j ' ) 

WRITE (6, 141) XM6 ,W6 _ _ 

1 4 1 FORMAT (T1 0 , ' XH6 ' ,719', ' = ' ,E 1 2 . 4,742,' MACH NUMBER , REGION 6 ' . /, 71 0 , 
*'W6' ,T19, '=' , E 12. 4, 742, 'JET WIDTH, REGION 6' ,TR2, ' (INCHES) ') 

WRITE (6,142) PT77DT77TT7 

142 FORMAT (TIP , 'PT7 ' ,719 ,'=' ,E12 . 4, T42 , 'TOTAL PRESSURE , REGION 7' ,742, 
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Figure B.3. — Continued 

A 


* • ,T82\ ' (PSIA) ' ,/,TlO, ' TT1 ' ;T19 ,Et2. 4/,T42, ’TOTAL TEMPERATURE , RE 

♦GIO N j 1 ,T82. ' (DEG.RANKIliE) ' -,r ' - = ■ -- . / ; 

WRITE (6,, 123) TT3 ,WA,WII 7" 

123 FORMAT (Tl 0 , 1 TT3 1 , Tl 9 , '=',E12.4,T42,' TOTAL TEMPERATURE REGION 3', 
*T82 , ' (DEG. RANKIHE ) 1 , /, T1 0 , ' WA *• , Tl 9 , ' = ' , E 12 . 4 , T4 2 /' MOLECULAR WEIGH 
*T ,A 1 R ' , T8-2 x ' ( 1/M OLE) ' , /,T1 0 , ' WH 1 ■, Tl 9 , « ,e 1 2 . 4 , T4'2 ; ' MOLECULAR W E IQ 

* IiT , HYDROGEN ' , T 8 2, ' (1/MOLE) ') / 

WR ITE ( G_, 1 24) Cl ,W 4,XM5 / 

12 4 FORMAT(Ti’0, '*Cl ' ,T19., ' = ' ,E12.4,T42, 'ORIFICE DISCHARGE COEFFICIENT' 

*J_ \_ _i «. '/■■■" • ■ 

♦,T10, 'W4' ,T19\' = ' ,E1 2 . 4 ,T42 , 'DIAMETER OF ORPFICE ' f T82 ^llNCHES)^' 
*/,TlO, 'XM5' ,Tl9v ' = ’ ,E12.4,T42, 'MACH NO. IN /REGION 5') 

IF (DEEP. LT. 0.0) "GO TO 202 / 

WRITE (6,125) PT1 , E)Tl ,T T1 / 

1 2 5 FORMAT ( // , T 1 0 , ' OUTPUT ',"// ,T10, 'PTl ' ,T19 7*=" r , E 12.4, T4T, 

* \ / 'TOTAL PRESS 

♦URE, REGION 1 ' ,T82, ' (PSIA) ' ,/,TlO, ' DT 1 '/, T 1 9 , ' = * , EX 2 . 4 T4 2 ~ 'TOTAL DE' 
♦NSITY , REGION 1 ' ,T82 , ' (LBM/CU . FT .) ' ,/,TlO ,.'TT1 < ,T19, ' = ' ,E12.4,T42, 
♦'TOTAL TEMPERATURE , REGION 1',T82,' (DEG. RANKINE) ') 

WRITE (6,126) ATI ,PT2 , DT2 / 

1 2 6 ' FORI LAT ( T 10,' AT 1 ' , T 1 9 , ' = ' \E 1 2 . 4 , T 4 2/ ' TOTAL SPEED OF SOUND , REG I ON T' ~ 

* , T82 , ' (FT. /SEC.) ' , /,T1 0 , * P V T2 ' , Tl 9/' = ' ,E1 2 . 4 ,T42 , ' TOTAL PRESSURE, RE 

♦Gl'Oli 2',T82, 1 (PSIA) ' ,7 , Tl 0\' DT2 ' /T 1 9 , ' = ' , E 1 2 . 4 ,T42 , ' TOTAL DENSITY, 
♦REGION 2 ' 7782., ( LBM/CU. FT . V ) J 

WRITE (6, 127) TT2 , AT2 ,P2 \“ / 

127 FORMAvT (T 1 0 . ' TT2 ' . T 1 9 , ' = ' , E 1 2 T 4 2 , ' TOTAL TEMPERATURE REG ION 2[ f 


I 


*T82 , ' (DEG . RANKIHE) ',/,TlO, t AT"2 ' ,T1 R , ' = ' ,E12.4,T42, 'TOTAL SPEED OF 

* SOUND , REGION _2 ' , T82, ' (FT./SEc\) ' ,/,TlO, 'P2' ,T19, ' = ' , E1 2 . 4 f T42 t 

* 'PRESSURE REGION 2' T82, * (PSIA)" 

WRITE (6, 128) D2,T2,A2 f 

1Z8 FORriATTT'T07'n^7%“TTT7^=^71TJ2T4 _ /T 4 2 ^DENSITY', REGION 2 T 7TR27 ' ( LV,K/CIT. 
♦FT.) ' ,/,TlO , ' T2 ' ,T19 , ' = ' ,<E 1 2 . 4 ,T 42 ^TEMPERATURE , REGION 2' ,T82, 

* ' (DEG . RANKINE ) ' , / , T 1 0 , 'JA2 ' , T 1 9 , • = ' >E 1 2 . 4 , T4 2 , ' SPEED OF SOUND , REG I 

♦ON 2' ,T82. ' (FT. /S EC.) ' )/ \ 

WRITE ( 6 , 1 2 9 j V2 , PT3 , DT6 ' V 

1 29 FORMAT (TIP, 'V2' ,T 19, ' = ' ,E12 . 4 . T42 , 'VELOCITY . REGION 2 ' ,T 82 , ' [FTj, /?. 
♦EC. j ' , / , T 1 0 , 1 PT3 ' , T 1 9 , ' = ' ,E12.4,T42, 'TOTAL PRESSURE REGION 3' , _ T82, 

* ' (PSIA) ' ,/,TlO, ' DT3 1 ,T19 , ' = ' ,E1 2 . 4 ,T42 ,\ TOTAL DEN SITY , RE GION 3' , 

* T82 , ' (LBM/CU .'FT. ) '/) ' ' 

WRITE ( 6 1 30 ) TT3 , AT- 3 , PT4 


130 FORMAT(T10, ' TT 3 ' f/ Tl9, ' = ' ,E12.4,T42, 'TOTAL 
♦T 82 ’ (DEG . RANKIHE ) ' ,/ ,TlO, ' AT 3 ' ,T 1 9, '=' , E 1 


k EMPERATURE , REGION 3 ' , 

. 4, T 4 2 ± 'TOTAL SPEED OF 


* SOUND, REGION 3/,T82, ' (FT./SEC.) ' ,/,T10,- ' PT4>^ , Tl 9,^ = T , El 2 . 4 ,T42 , 

* '.TOTAL PRESSURE. REGION 47, T 82, ' (PSI A) ' ) 


VJRITE (6,131) DT4 , TT4 , AT4 

131 FORMAT (T10,/PT4' ,T19, '=' ,E 12.4,T42 , 'T OTAL DENSI TY, REGION 4' ,T82. 

* ' ( LBM/CU. F^T. )' , / , T 1 0 , ' TT 4 ' , T 1 9 , ' = ' ,E12.4,T42, 'TOTAL TEMPERATURE , R 
♦EG ION 4 7^7,8 2, ' (DE G. RANKINE) ' ,/,T10 , ' AT4 1 ,T 1 9 , '=^, E 12 . 4,T4 2 , ' TOTAL 

♦ SPEED OF "’SOUND , REGION 4' ,T82, ’ (FT./SEC.) ’)” 

WRITE (6. 132) P4.D4.T4 
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Figure B. 3. — Concluded ! 

* ' (PSIA) ' ,/,Tl 0 , 'DT7' ,T19, '=' ,E12.4,T42, ' TOTAL - ' DENSITY , REGION ' 7', 

/CII. F2M. ' . /.TIP. 'TT7' .T19. '=' .E12 . 4 ,T42 . 'TOTAL TEMPERAT 

*URE, REGION 7 ' ,T82 , ' (DEG. RANKINE) ' ) ' •" 

W RI TE (6_». 14 3). AT7 ,P7 ,D7 . 

143 FORMAT (T1 0 , 1 AT7 ' ,T1 9 , '=' ,E12.4,T42, 'TOTAL SPEED OF SOUND , REGION 7' 
*^82.,, ' (FT. /SEC. ) ' . / T 1 0 . * P7 1 .TIP . ' = ' .E12.4.T42. ' PRESSURE . R EGION 7 ' . 

* T82 , ' (PSIA) ' , /T1 0 , ' D7 ' ,T .1.9 , /=-' ,E1 2.:4 ; ,T42 , '-DENS IT Y, REGION 7 ' , 

_T 82 , l (LBM . /CU. FTi) ' ) • Li -• 

WRITE (6, 144) T7,A7,V7 

1 4 4 FORMAT (TlO^ 'T7.1 .T IP . 1 = ' ,E12. 4.T42 . ' TEMPERATURE , REGION 7 ' ,T 82, __ 

* • (DEg". RANKINE) ' , / , T 1 0 , ' A7 ' „T1P , ' = ' , E 1 2 , 4 , T4.2 , SPEED OF SOUND, REGl' 
♦ON 7 » , T8 2 . ' (FT . / SEC . ) ./.TIP . 1 V7 ' . TIP . ' ± ' ,E T2 . 4 . T4 2 ; 'VEL OC ITY, RE G I 
♦ON 7',T82, ' (FT./SEC.) ') 

WRITE (6 . 1 _4 5 ) XM 7.W7 ] ■ 

145 FORMAT (TIP, ' XM7 * , T 1 P , '= ' ,E12 . 4 ,T42 ,,'MACH NUMBER, REGION 7 ' ,/,T 10, 

* ]yn ' , T 1_P J = ' , E 12.4, T 42, *JET WIDTH, REGION 7'. ,T82, ' (INCHES )' ' ) ’ 7. 

WRITE (6 ,146) F4,F4P,F5 

1 4 G _ FORMAT (T 10 , ' F 4 ' ,TlP , ' = 1 , E 1 2 . 4 , T42 , 'MASS FLOW, REGION 4',T 82, ' (LBM . / 
♦SEC, f ' ,/,'tTO , 'F4P ' ,T~1 9,, ,E1 2 . 4 , T 42 EFFECTIVE MASS FLOW REGION'" V 

* ' ,T32, ' LBM . /SEC . )' . / ,T1 0 . *F5 * .TIP. &&E 12 . 4 . T42 .'MASS -FLOW , REGI ON ' 

* 5 ’ , T82 , ' (LBM . /SEC . ) 1 ) 

_ WRITE ( 6 ,J_4 7 ) F6 # F7 . QR ; . 

147 FORIIAT (T 1 0 , ' F6 ' ,Tl9 , ' = ' ,El 2 . 4 ,T42 , -'MASS FLOW, REGION , 6 V,T82 , ' (LBM. /" 
♦SEC. ) ' ,/,T 10, ' F7 '_*T1 9 . ' - * .El 2. 4. T 4 2. 'MASS FLOW. REGION 7 ' ,.T 82,' (LB M 
* . /SEC. ) ' , / , T 1 6 , ' QR ' ,T19, '=' ,E12. 4,T«2 , ' DYNAMIC PRESSURE RATIO ' ) 

_ WRITE (6. IIP) S5 QS4P 

IIP FORMAT (T10, 'S50S4P' ,T19, ' = ' , E 1 2 . 4.,T4:2 , 'AREA RATIO REGION 5 TO REGl" 

♦on . 4p '_) . • ; X ;■ r >./ ; ; . . 

GO TO 22 

_ 20.2 WRITE (6 L 20 3 )__ _ 

203 FOlVlAT MOX, 'DEEP LESS THAN ■ 0.0 J) ,. ; , v y ' 

1 STOP . 111 ; . ' . 

END 

_/e.ATA _ 

1.4 4.03 300.3 

_ 294. G 517. __ 509. 

29.0 2.016 0.5 

.040 

4.6 


.999 

END 
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Figure B.4. Sample Output from Barrel Shock Computer Program 


INPUT 





G 


0.1400E 01 

SPECIFIC HEAT RATIO 

. ■ - 

XM2 

= 

9 . 40 30)1 01 

MACH HO. Ill REGION 2 , E P.F.F.STPT1 AM 

' 

PT 1 


0.300 2E 03 

TOTAL PRESSURE , REGION 1 

(r»SIA) 

I:r:i -1 

i .4. ..1 

= 

0.2046]-; 0 2 

TOTAL PRESSURE REGIOU2 

;(r»cTA) 

TT 1 

=2 

0. 5170E 03 

TOTAL 'TEMPERATURE , REGION 1 

t>a?*KI?!F) 

.J. .L J 

= 

0.50001; 03 

1 TOTAL TEMPERATURE REGION 3 . 

1 (r)r»r* r nT'.TTjf yrTT7^ 

s ; A 

= 

o . 200 or; o? 

I !- ’.OLE Cl -LAP WEIGHT , A 1 R 

( 1 /M.OT/E ) 

t n r 

= 

0.2016 r; oi 

MOLECULAR WEIGHT , HYDROGEN 

( 1 /MOLE) 

Cl 

=5 

0.5000E no; 

ORIFICE DISCHARGE COEFFICIEET 


174 

= 

o.40ooi:-oi 
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Figure B.3. — Concluded 
• _• 
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* TR2 f ! (LBH . /CU .FT.) 1 . /.TIP. ' TT7 ' .T19. ' = 1 . E 1 2 . 4 , TJi 2 , ' TOTA L TEMPERAT 

♦URE/EEGION 7' ,T82, ' (DEG. RANKINE) ') ‘ *' / 

WR ITE (6 f 14 3 ) A T 7.P7.D7 ■ £__ 

143 FORMATs(TlO , 'AT7' ,T19, ' = ' ,E12.4,T42, 'TOTAL SPEED OF SOUND, REGION 7' 

* .T82 . 1 (FT. /SEC. ) 1 ,/TlO , 'P7' .TIP . ' = ' ,E1 2 . 4 ,T4"2 . ’PRESSURE . R EGION 7 ' . 

* T82, ' (PSIA) ' ,/TlO, 'D7',T1.9, ' = ' ,E-1 2'.4.;,T42 , ^DENSITY, REGION 7', 

_* £82 (LBM . /CU . FT.) ' ) ' ■ / f ' _ " / : - _______ L 


WRITE (6, 14 s 4) T7 ,A7 ,V7 / 

144 FORMAT (T1oAt7| ,T19 . ' = ' ,E12. 4.T42. 'TEMPERATURE, REG ION 7* ,T82, 

"'■*T(6eg 7' RANKINE) ' ,/,TlO , '-A7' ,T19 , ' = ' ,E12’. 4,T4.2 , ' SPEED OF SOUND, RE'gT 

♦ ON V ,T82 (FT. /SEC.) 1 ./.TIP. 'V7* .Tt9/ , = t .El 2 . 4 ,T42 j ' VEL OC ITY, RE GI 

♦ON 7 ' ,T82 , ' (FT 7/ SEC. ) * ) / 

WRITE (6.1 _4 5 ) XM 7\.W7 / : 

1 4 5 FORTIAT (T1 0 , 1 XM7 ' /T 1 9 , ' . = ' , E 1 2 . 4 , T 4 2 * MACH NUMBER , REG ION 7 1 , / , T 1 0 , 

♦ ' W7 ’ 


MAT ( T 1 0 , ' XM7 * , N T19, ' = ' ,E 1 2 . 4 ,T424'!MACH NUMBER, REGION 7',/,T 

' . T 1 9 . 1 = ' . E 1 2 . 4\T 4 2 , ' JET WIDTH./REGION 7* ,T8 2 , * ( INCHES )' * ) . 

WRITE (6, 146) F4,F4PVF5 / ‘ 
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1 sto p J_ 

END 
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Figure'B.4.'— Concluded 
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APPENDIX C 


MOLECULAR VISCOSITY FOR H 2 /N 2 /0 2 MIXTURES 


In modeling of turbulent transport coefficient distributions, it is necessary to know the 
molecular viscosity at the wall. This parameter, p m j x , * s usec * * n von Priest’s damping function, 
and is also the limiting value for the transport coefficient of momentum at the wall where e = 0, 
i.e., 

Meff = e + Fmix ( QA) 

where e is the turbulent eddy viscosity. 

For flow containing mixtures of gaseous hydrogen, oxygen, and nitrogen, the molecular 
viscosity may be calculated as follows (ref. 33): 

N 

^mix - ^ x i^i^i (C.2) 

1=1 


fi 

where 

<Ajj 



(C.3) 


(C.4) 


Wj = Molecular weight of species i 
Pj = Viscosity of species i 
Xj = Mole fraction of species i 

N = Total number of species considered 

A computer program was written (fig. C. 1) to solve equation (C.2) for M ni j x , and calculat- 
tions were performed to check the subroutine and determine the dependence of M m j x on mass or 
mole fraction of hydrogen. The program was debugged against hand-calculated and redundant check 
cases. Figure C.2 shows the computed dependence of P m j x of a hydrogen/air mixture. The high 
degree of nonlinearity between M m j x and xj-j^ is consistent with comments in reference 33. 
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H. 000 1 
L. 0002 
L . 0003 
L.0004 
L. 0005 
L . 0 00 6 
L.0007 
L . 0 0 0 3 
L . 0009 
L.0010 
L.0011 
L.0012 
L. 00 13 
L.0014 
L.0015 
L.001G 
L.0017 
L.0010 
L.0019 
L. 0020 
L.0021 
L . 0022 
L. G 023 
L.0024 
L. 0025 
L.002G 
.Li • 0 0 2 7 
L. 0020 
L . 0029 
L. 0030 
L.0031 
L . 00 32 
L. 0033 


SUERO UT I HE XMUMI X ( X , XI IUBAR ) 

C CALCULATE VISCOSITY IN II2/02/N2 GAS MIXTURES 
C ENTER ROUTINE WITH MOLE FRACTIONS X { I) . . . 1=1 (K2 ) , 2 (02 ) , 3 (N2 ) 
C LEAVE WITH VISCOSITY XMUBAR (LBM/FT-SEC) 

C 

C 

C SEE PG. 24 or BIRD, STEWART AND LICHTFOOT 

C FOR 112 , 02 , M2 MIXTURES 1=1, 2, 3 RESPECTIVELY 

DIMENSION 1/(3) ,XMU(3) ,X(3) , PHI ( 1 0 , 1 0 ) ,F(10) 

DATA W/2 . 0 1 6 , 32 . , 2 8 . 0 1 G/ 

DATA XMU/. 00035, .0189, .01G57/ 

C XMU IN CENTIPOISE 

DATA I ONCE/ 0/ 

IF (IONCE) 4,4,5 

4 IONCE= 1 

mi (1 ,2) =1 . 05G4 
1*111(1 , 3) =1 .919 
PEI (2, 1 ) =.2G472 
F 111 (2,3)=. 99347 
PHI (3, 1) =.27403 
PHI ( 3 ,2) = . 9993G 
PHI (1,1) = 1 . 

PEI (2,2) = 1. 

PHI (3,3) = 1. 

5 CONTINUE 
DO 1 1=1,3 

1 F ( I ) = X(1 ) *PIII (1, 1 )+X(2) Till (I, 2) +X (3) +PHI (I, 3) 

XMUBAR = 0. 

DO 2 1=1,3 

2 XMUBAR '= XMULAR+X (I) *XMU (I) /F (I) 

XMUBAR = XMUBAR * 6. 7 2 E- 4 

RETURN 

END 


Figure C. 1 . Listing of Mixture Viscosity Computer Program 
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APPENDIX G 

MOLECULAR VISCOSITY FOR H 2 /N 2 /0 2 MIXTURES * 

/ 



In modelin^of turbulent transport coefficient distributions, it is' necessary to know the 
molecular viscosity autfhe wall. This parameter, M m j x , is used in von^Driest’s damping function, 
and is also the limiting value for the transport coefficient of momentum at the wall where e = 0, 


i.e., 


^eff e+ ^mix 



(C.l) 


where e is the turbulent eddy viscosity. 

\ 

For flow containing mixturesv^f gaseous hydrogen, oxygen, and nitrogen, the molecular 

calculated as follows 6ref. 33): / 

/ 


viscosity may be 


mix 


where 


0 ; 



N 

2 Xj/ij/fj 


N 

2 0ij xj 



(C.2) 


(C.3) 


/Mj\ /Z / w\ 1/4 -i 2 

W vv\ 


W: = Molecular weight of species i 


= 

x i = 
N = 


Viscosity o^pecies i 
Mole fraction of species i 
Total number of species considered 
/ 



(C.4) 


A computer program was written (fig. C. 1) to solve equation (C.2) Tor M m j x , and calculat- 
tions were performed to check the subroutine and determine the dependence, of M m j x on mass or 
mole fraction of hydrogen. The program was debugged against hand-calculated and redundant check 
cases. Figure C.2 shows the computed dependence of M m j x of a hydrogen/air mixture. The high 
degree of nonlinearity between /i m j x and x^j is consistent with comments in reference 33. 

\ 
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